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Abstract 

This  work  investigates  the  aggregate  performance  of  satellite  receiver  formations 
functioning  as  orbiting  interferometers  as  compared  to  filled  apertures  of  similar  geome¬ 
tries.  The  resulting  models  facilitate  selecting  initial  conditions  for  formations  such  that 
their  control-free  dynamics  yield  interferometry  performance  with  near-minimal  errors  as 
compared  to  the  filled  apertures.  The  solution  method  draws  on  the  dynamic  models  of 
an  orbiting  planar  satellite  formation  to  define  the  size  and  shape  of  a  reference  aperture 
and  to  define  the  degrees  of  freedom  for  the  formation  members.  The  modelled  paths 
of  formation  elements  yield  predictable  geometries  at  any  time  for  which  the  aggregate 
performance  of  the  array  of  discrete  receivers  may  be  calculated.  The  objective  of  the 
optimization  process  is  to  minimize  the  time-averaged  square  of  the  difference  between 
the  filled  aperture’s  intensity  map  and  that  generated  by  the  discrete  receiver  array.  This 
yields  a  formation  whose  configuration  offers  near-minimum  errors  for  imaging  processes 
beginning  at  any  arbitrary  start  time.  The  problem  as  posed  is  non-convex,  and  requires 
implementation  of  a  global  search  method.  Genetic  algorithms  are  used.  The  genetic  al¬ 
gorithms  construct  populations  of  generic  coefficients,  adaptable  to  the  number  of  degrees 
of  freedom  for  each  proposed  formation.  The  fitness  assigned  to  members  of  each  popu¬ 
lation  is  the  average  aggregate  error  its  formation  generates  over  the  simulation  duration. 
The  solution  method  includes  a  new  analytic  solution  for  the  intensity  map  of  an  elliptical 
aperture  and  a  technique  for  generalizing  this  solution  to  include  the  effects  of  non-ideal 
viewing  geometries. 
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MEAN  PERFORMANCE  OPTIMIZATION  OF  AN  ORBITING  DISTRIBUTED 


APERTURE  BY  WARPED  APERTURE  IMAGE  PLANE  COMPARISONS 

I.  Introduction 

1.1  Overview 

The  notion  of  creating  satellite  formations  for  distributed  missions  is  by  no  means 
new,  but  has  received  renewed  interest  in  the  past  few  years.  Distributing  a  satellite 
mission  over  a  formation  of  smaller  vehicles  offers  the  potential  advantages  of  very  large 
(although  sparse)  effective  apertures,  modular  maintainability  and  upgradability,  graceful 
degradation,  and  greatly  reduced  life  cycle  costs. 

Distributed  satellite  systems  (DSSs)  are  among  the  most  challenging  conceptual  de¬ 
sign  problems,  since  modelling  such  systems  requires  nonlinear,  high  order,  often  coupled 
mathematics  [55,  54,  18].  The  Air  Force  and  NASA  spurred  feasibility  research  in  the 
creation  of  distributed  satellite  missions  through  the  TechSat  21  research  initiative.  The 
TechSat  21  program  encouraged  investigation  into  not  only  satellite  formation  dynam¬ 
ics,  but  also  micro-satellite  and  micro-propulsion  design,  distributed  mission  architecture, 
sparse  aperture  sensing,  collaborative  behavior,  and  micro-nano-technology  [41].  These  in¬ 
vestigations  resulted  in  tremendous  advancements  in  the  science  of  modelling  the  dynamics 
of  satellite  formations  and  optimizing  control  use  for  their  maintenance  and  reconfigura¬ 
tion.  Other  work  has  addressed  the  potential  uses  of  formations  for  communications,  radar, 
astronomy  [6,  17],  and  moving  target  identification  [38,  39].  Many  of  these  works  addressed 
optimization  approaches  for  formations  used  for  interferometry  missions  [9,  20,  23,  7]. 

This  interdisciplinary  work  draws  from  the  fields  of  astrodynamics,  physics,  electrical 
engineering,  and  mathematics.  It  melds  the  disparate  sciences  of  orbital  mechanics,  optics, 
and  optimization  to  build  a  solution  methodology  for  shaping  satellite  formations.  The 
basic  nature  of  formations  makes  them  attractive  platforms  for  receiver  arrays  performing 
interferometry.  Shaping  a  formation  to  optimize  its  performance  is  necessarily  both  an 
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optical  problem  and  an  orbital  mechanics  problem.  Orbital  mechanics  restricts  the  design 
choices  for  formation  shapes  and  dictates  how  receivers  move  within  that  formation.  The 
science  of  optics  defines  both  the  performance  of  a  filled  aperture  of  size  comparable  to  the 
formation  and  the  aggregate  performance  of  the  formation’s  discrete  receivers.  Building 
these  elements  into  the  fitness  model  offers  a  methodical  approach  to  minimizing  the  dif¬ 
ference  between  the  outputs  of  these  two  apertures  by  adjusting  the  defining  parameters 
of  the  formation’s  members. 

Optimal  geometries  are  a  key  element  in  the  employment  of  distributed  receivers  for 
astronomical  observations  [42,  28],  but  are  not  limited  to  stellar  targets.  Optimizing  re¬ 
ceiver  arrays  for  point  sources  is  a  logical  approach  for  maximizing  performance  in  imaging 
very  small  targets.  Literally  dozens  of  metrics  for  performing  formation  optimization  have 
been  designed.  Some  minimize  control  usage  or  reconfiguration  time.  The  thesis  and  series 
of  papers  by  Sedwick,  Miller,  and  Kong  at  MIT  [21]  optimized  instantaneous  performance 
of  a  formation  by  comparison  to  a  filled  aperture.  This  work  generalizes  this  approach 
to  minimize  errors  over  an  arbitrary  imaging  time.  This  approach  requires  taking  into 
account  the  loci  swept  out  by  the  formation  receivers  by  shaping  the  reference  aperture  to 
match  these  loci.  Therefore  the  control- free  dynamics  of  a  formation  are  used  to  define  an 
appropriately  shaped  aperture,  and  candidate  formations  are  compared  to  this  baseline. 
The  intensity  map  generated  by  a  simple  circular  aperture  is  a  well-known  solution  and  is 
the  reference  nearly  always  used  in  previous  works.  Such  a  reference  is  inefficient  since  it 
inherently  either  includes  receiver  paths  not  unreachable  due  to  the  formation  dynamics  or 
ignores  beneficial  paths  that  are  reachable.  This  work  includes  a  new  solution  method  for 
mapping  physical  distortions  of  a  circular  aperture  from  the  physical  or  spatial  frequency 
domain  through  the  Point  Spread  Function  into  the  spatial  domain  of  the  image  intensity 
map.  Developed  first  for  an  ideal  viewing  geometry,  this  solution  method  is  then  general¬ 
ized  to  include  the  apparent  aperture  distortions  generated  by  realistic,  non-ideal  viewing 
geometries. 

The  optimization  task  is  therefore  one  of  selecting  the  defining  parameters  of  a  for¬ 
mation  such  that  the  the  average  difference  between  a  reference  aperture’s  intensity  map 
and  that  generated  by  the  aggregate  performance  of  the  discrete  receivers  is  minimized. 
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This  turns  out  to  be  a  non-convex  problem.  The  lack  of  convexity  is  overcome  by  making 
a  few  simplifying  assumptions  and  exploiting  the  powerful  global  search  capabilities  of 
genetic  algorithms.  The  defining  parameters  of  formation  dynamics  limit  the  number  of 
degrees  of  freedom  in  a  series  of  increasingly  complex  problems,  each  optimized  by  genetic 
search  routines.  As  the  simplifying  assumptions  are  relaxed,  the  search  procedure  selects 
a  growing  set  of  generic  coefficients  each  adaptable  to  the  varying  dimensionality  of  the 
increasing  degrees  of  freedom.  With  greater  numbers  of  receivers  in  the  formation  and 
increasing  flexibility  in  their  initial  conditions,  the  discrete  array’s  ability  to  synthesize  the 
performance  of  a  filled  aperture  grows  rapidly.  As  the  formation  size  grows,  the  aggregate 
error  dwindles  to  near-zero  indicating  that  the  shaped  aperture  reference  is  indeed  a  valid 
ideal,  theoretically  reachable  by  an  infinite  number  of  receivers. 

1.2  Research  Contributions 

Though  this  work  relies  heavily  on  the  pioneering  efforts  of  many  experts  in  the  fields 
of  optics  [13],  interferometry  [46,  35],  orbital  mechanics  [50],  and  formation  dynamics  [56, 
4] ,  it  offers  significant  advancement  in  several  aspects  of  both  formal  problem  definition  and 
adaptation  of  global  search  algorithms  for  its  solution.  The  purpose  of  this  research  is  to 
encapsulate  the  relevant  dynamics  of  satellite  formations  and  the  physics  of  interferometry, 
specifically  as  it  relates  to  receiver  arrays  both  liberated  and  restricted  by  placement  in 
orbit.  By  applying  the  science  of  optics  to  both  an  aggregate  aperture  and  a  filled  aperture, 
a  metric  is  defined  as  a  comparison  of  the  two,  then  expanded  to  include  time  variations. 
Employing  this  as  the  measure  of  merit  for  simplified  formations  gives  rise  to  a  non-convex 
problem.  The  global  search  techniques  of  genetic  algorithms  are  well  suited  to  non-convex 
problems  and  provide  acceptable  numerical  solutions.  The  specific  contributions  stemming 
from  this  problem  definition  and  solution  are  as  follows: 

i.  The  model  built  in  this  research  characterizes  the  process  for  building  an  intensity 
map  generated  by  an  aggregate  aperture  of  discrete  receivers,  Ip(£,  t),  and  formal¬ 

izes  the  mathematics  for  building  a  complete  set  of  baselines  for  a  generalized  set  of 
asymmetrically  placed  receivers. 


1-3 


ii.  The  second-order  instantaneous  metric  of  previous  work  [20]  is  extended  to  include 
variations  over  a  finite  dwell  time.  This  metric  holds  the  promise  of  yielding  forma¬ 
tions  whose  control-free  motion  yields  near-minimal  imaging  errors  for  any  arbitrary 
start  time. 

iii.  The  mathematically  simple  circular  reference  aperture  is  discarded  in  favor  of  more 
complex  shapes  which  cover  all  available  receiver  loci  and  include  none  that  are 
unreachable  due  to  the  formation  dynamics. 

iv.  Both  a  numerical  and  an  analytic  solution  are  offered  for  the  intensity  map,  I0(£,  rj), 
of  a  filled  elliptical  aperture  of  arbitrary  eccentricity.  An  aperture  with  this  shape 
covers  all  relative  loci  of  formation  elements  in  control-free  motion  (as  viewed  from 
directly  above) .  The  intensity  map  offered  by  the  analytic  solution  is  shown  to  match 
the  less  efficient  but  mathematically  simpler  numerically  derived  map. 

v.  The  analytic  solution  for  an  elliptical  aperture  is  extended  to  include  the  distorting 
effects  of  general  viewing  geometries.  These  geometries  are  determined  by  the  Kep- 
lerian  elements  defining  the  reference  orbit  and  by  the  inertial  locations  of  targets. 

vi.  The  generalized  aperture  solution  demonstrates  a  new  technique  for  generating  an 
analytic  intensity  map  for  any  (circular  or  elliptical)  aperture  distorted  as  viewed 
from  general  locations.  The  solution  method  may  readily  be  adapted  for  use  with 
any  orbital  reference  or  target  mapping  convention. 

vii.  Global  search  techniques  are  adapted  to  generate  generic  coefficients  which  select 
the  defining  parameters  of  simplified  formations.  The  generic  coefficient  approach 
allows  for  the  selection  of  larger  sets  of  defining  parameters  of  each  formation  element 
without  significant  code  rewriting.  The  varying  dimensionality  of  these  parameters 
is  likewise  transparent  to  the  selection  routine. 

1.3  Dissertation  Outline 

This  dissertation  consists  of  8  chapters,  including  this  introduction.  Chapter  2  re¬ 
views  the  extensive  literature  on  modelling  the  dynamics  of  satellite  formations  and  the 
limitations  of  several  leading  model  solutions.  Each  model,  despite  extensive  differences 
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in  their  formulations,  is  shown  to  limit  the  formation  designer’s  choices  to  a  maximum  of 
four  degrees  of  freedom  per  element. 

Chapter  3  summarizes  the  basic  science  of  interferometry  from  the  nature  of  a  single 
pair  of  interfered  wave  forms  out  to  the  synthesis  of  an  entire  aperture’s  intensity  map. 
Chapter  4  reviews  the  extensive  literature  on  moving  interferometry-based  receivers  out 
to  orbital  platforms  and  the  assumptions  and  limitations  inherent  in  the  various  system 
descriptions.  The  analytic  solutions  for  aggregate  aperture  performance  are  re-derived 
for  symmetric  arrays  to  provide  background  for  the  more  general  formalized  treatment 
of  asymmetric  arrays.  These  are  incorporated  into  the  development  of  a  tinre-averaged 
comparison  metric.  It  formalizes  the  aggregate  aperture  construction,  the  selection  of  an 
appropriately  shaped  reference  aperture,  and  determines  a  conservative  size  for  the  image 
plane  over  which  the  comparison  is  performed. 

Chapter  5  takes  the  basic  technique  of  integrating  a  point  spread  function,  PSF, 
over  an  aperture  to  generate  its  intensity  map  and  applies  it  to  apertures  with  simple 
geometries.  Beginning  with  a  simple  one-dinrensional  aperture  (a  slit),  it  is  then  applied 
to  simple  two-dimensional  shapes.  The  elliptical  aperture  derivation  at  the  end  illustrates 
the  basic  technique  for  warping  a  circular  intensity  map  in  one  direction,  a  technique 
generalized  in  chapter  7. 

The  nature  of  the  optimization  problems  based  on  the  metrics  of  chapter  4  is  ad¬ 
dressed  in  chapter  6.  It  begins  with  the  simplest  possible  one-dinrensional  formation  to 
illustrate  the  convexity  of  the  problem  and  gradually  expands  the  trial  formations  to  more 
realistic  relative  geometries  and  more  freedom  for  receiver  loci  selection.  Chapter  7  finally 
adds  in  the  apparent  distortions  caused  by  viewing  formations  from  non-ideal  directions. 
These  distortions  affect  both  the  receiver  loci  and  the  filled  aperture.  Four  example  stars 
are  used  to  demonstrate  the  nature  of  each  step  in  the  process  of  distorting  the  apparent 
loci  of  receivers  about  a  reference  orbit  and  mapping  the  distortions’  effects  to  the  im¬ 
age  plane.  The  elliptical  aperture  solution  from  chapter  5  is  expanded  to  likewise  include 
these  distortions  to  a  filled  aperture.  This  solution  is  incorporated  as  the  reference  map, 
I0(£,  77)*,  in  the  metrics  of  chapter  4  and  several  example  formations  are  optimized  using 
the  techniques  of  chapter  6.  The  optimized  formations  included  there  and  in  Appendix  F 
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illustrate  the  aggregate  error  approaching  zero  for  formations  with  increasing  numbers  of 
receivers.  The  remaining  appendices  describe  known  solutions  for  formation  dynamics  and 
interferometer  performance  all  adapted  to  the  standardized  notation  of  this  research. 
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II.  Satellite  Formation  Dynamics 

Describing  the  dynamics  of  a  satellite  formation  relative  to  a  localized  frame  is  a  surpris¬ 
ingly  difficult  problem.  The  models  used  to  describe  the  motion  of  formation  elements 
determine  both  the  available  degrees  of  freedom  for  choosing  their  configurations  and  the 
relative  loci  that  receivers  trace  out  over  time.  This  chapter  begins  with  the  earliest  treat¬ 
ment  of  this  problem  by  Clohessy  and  Wiltshire,  as  expanded  by  Hill,  and  improved  by 
Schaub  and  Alfriend  to  reveal  the  constraints  on  element  paths.  The  final  section  describes 
the  latest  modelling  improvements  by  Wiesel  that  vastly  improve  the  fidelity  of  element 
propagation,  but  still  exhibit  the  four  degree  of  freedom  limitation  of  earlier  models. 

2.1  Clohessy- Wiltshire  Dynamics 

The  dynamics  of  satellite  formations  in  a  local  coordinate  system  had  its  first  useful 
treatment  by  W.  H.  Clohessy  and  R.  S.  Wiltshire  [5]  back  in  1960,  as  they  sought  a 
simple  solution  to  be  used  for  short  time-frame  orbital  rendezvous.  Clohessy  and  Wiltshire 
simplified  the  problem  by  linearizing  the  orbital  dynamics  about  an  unperturbed  circular 
orbit.  This  circular  orbit  provided  the  origin  for  a  local  coordinate  system,  with  axes 
directed  along  the  radial,  velocity1,  and  orbit  plane  orthogonal  directions.  The  first-order 
linearization  of  the  orbital  dynamics  in  this  reference  frame  yields  a  constant  coefficient 
system.  The  resulting  solution,  see  Appendix  A, 


x  =  2coy  +  3uj2x 

y  =  —2  tax  (2.1) 

z  =  —u2z 


hereafter  referred  to  as  the  C-W  solution,  has  provided  at  least  the  initial  setup  for  most 
investigations  into  formation  dynamics  and  control. 

The  basic  C-W  solution  is  an  adequate  model  for  short-term  applications  (hours), 
and  has  been  used  for  rendezvous  studies.  Guelman  and  Aleshin  [12]  used  the  basic  linear 


lrThe  original  paper  by  Clohessy  and  Wiltshire  actually  chose  the  antivelocity,  radial,  and  orbit  normal 
directions  for  the  local  coordinate  system.  The  more  conventional  choice  of  radial,  velocity,  and  orbit 
normal  directions  is  used  throughout  this  work. 
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dynamics  to  generate  rendezvous  maneuvers  that  minimized  fuel  usage  with  arbitrary 
constraints  on  thrust.  Unfortunately,  since  the  linearized  solution  does  not  incorporate 
even  the  J2©  perturbation,  it  lacks  sufficient  fidelity  for  long  term  control  law  synthesis. 
The  higher  order  effects  not  modelled  by  the  C-W  solution  incite  unacceptably  high  control 
usage  and  reduce  satellite  life  spans  to  impractical  durations. 

Schaub  and  Alfriend  [37]  extended  the  C-W  solution  to  include  first  order  oblateness 
effects  (effects  of  order  ,/|^  or  higher  are  ignored).  Their  solution  described  the  relative 
orbits  in  terms  of  the  Delaunay  orbital  elements  [50]  ( l ,  g,  h,  L,  G,  H ),  where  l  is  the  mean 
anomaly,  g  is  the  argument  of  perigee,  h  is  the  right  ascension  of  the  ascending  node,  and 
(L,  G,  H )  are  the  associated  generalized  momenta.  Since  the  first  order  J2©  effects  are 
not  a  function  of  the  mean  orbital  elements  (l,g,h),  these  are  unconstrained.  It  is  not 
possible  to  eliminate  the  drift  rates  of  each  orbit  element,  so  Schaub  and  Alfriend  sought 
in  their  formulation  to  keep  the  differences  in  mean  orbit  element  growths  to  zero,  thus 
eliminating  relative  secular  growth.  However  eliminating  relative  secular  rates  places  two 
constraints  on  the  choices  of  (l,g,  h ).  Choosing  any  of  the  first  three  orbital  elements  and 
any  one  of  the  momenta  completely  determines  the  other  two  momenta.  This  creates  a  four 
degree  of  freedom  system  rather  than  six  for  each  formation  element.  Eliminating  all  of  the 
drift  rates  actually  places  a  third  constraint  on  the  momenta,  completely  determining  the 
formation  from  the  angular  elements  and  severely  limits  the  choices  for  formation  design. 

The  orbital  period  constraints  on  each  element  are  more  readily  seen  in  Hill’s  param¬ 
eterization  [14]  of  the  C-W  solution.  Appendix  B  summarizes  the  description  of  formation 
orbits  in  terms  of  six  constants  (aH.  bH,  pH,  mH,  nH,  and  0H,),  of  which  aH  must  be  0 
for  each  element  to  eliminate  secular  growth  in  its  path.  The  parameter  bH  is  simply  an 
offset  of  the  entire  formation  along  the  velocity  vector  direction,  ey.  Since  introducing  a 
constant  offset  in  the  ey  direction  adds  nothing  to  interferometry  performance,  without 
loss  of  generality  this  parameter  is  assumed  0  for  each  element.  Again,  this  yields  four 
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degrees  of  freedom  per  formation  element. 


x(t)  =  pH  sin (cat  +  9h)  +  aH 

y(t)  =  2pHcos(ujt  +  0H)  -  ^faHt  +  bH  (2.2) 

z(t)  =  mHpH  sin(o;t  +  0H)  +  nH2pH  cos (cat  +  0H) 

Yeh  and  Sparks  [56]  analyzed  Hill’s  equations  to  develop  geometries  of  free-flying  satellite 
formations  without  the  use  of  controls.  Control-free  formations  limited  their  design  degrees 
of  freedom  to  four  per  satellite,  in  agreement  with  the  initial  formulations  of  Schaub  and 
Alfriend.  Inspection  of  equation  2.2  reveals  that  each  element  in  the  formation  traces  out 
an  elliptical  path  with  a  2:1  ratio  of  major  to  minor  axes. 

A  variety  of  optimal  control  schemes  have  been  developed  using  the  limited  fidelity 
solutions  discussed  above.  Yeh,  Nelson,  and  Sparks  [55]  developed  a  methodology  for 
generating  control  laws  using  a  sliding  mode  framework.  Their  work,  based  on  the  basic  C- 
W  solution,  incorporated  J2©  effects  by  increasing  the  controller  bandwidth.  The  thesis  by 
Irvin  [15]  investigated  minimum  fuel  reconfiguration  techniques  using  the  C-W  solution  and 
a  variety  of  feedback  design  techniques  (linear,  LQR,  SDRE,  sliding  mode).  He  found  that 
formation  reconfiguration  could  be  accomplished  for  minimal  fuel  usage  using  the  simplest 
of  linear  techniques.  Chichka  [4]  modelled  satellite  formations  to  remain  in  a  constant 
apparent  distribution  as  viewed  from  the  planet  surface ,  not  an  inertial  or  co-orbiting 
reference.  His  selection  of  formation  geometry  was  strongly  driven  by  viewing  angles  from 
a  surface  location,  but  did  demonstrate  the  feasibility  of  creating  the  formations.  DeCou  [9] 
used  the  basic  C-W  solution  to  develop  a  station-keeping  strategy  for  three  satellites  forced 
to  remain  in  the  (u  v)  plane  (see  Section  3.2  for  a  description  of  the  {u  v)  plane).  This  plane 
was  to  be  oriented  at  an  arbitrary  angle  to  the  orbital  plane,  specifically  for  performing 
interferometry.  The  scheme  assumed  a  constant  thrust  (via  ion  engines)  available  for  the 
control  system’s  use,  and  sought  to  minimize  propellant  usage. 

Sedwick,  Miller,  and  Kong  [40]  used  the  dimensional  equation  techniques  of  Bucking¬ 
ham  [3]  to  include  more  perturbation  effects  into  Hill’s  equations.  The  dimensional  equa¬ 
tions  cite  that  a  system  with  <7  characteristics  of  interest  with  dimensionality  77,  should  be 
describable  by  a  —  i]  dimensionless  parameters.  Looking  only  at  planar  arrays,  this  pro- 


2-3 


duced  two  descriptive  constants  (similar  to  the  four  degrees  of  freedom  from  Appendix  B, 
if  rnH  =  nH  =  0).  They  found  that  first  order  J20  effects  were  far  stronger  than  any  other 
perturbations,  even  out  past  geosynchronous  altitudes.  The  differential  J2©  effects  were 
comparable  with  solar  radiation  pressure,  and  the  two  effects  dominated  the  perturbations 
of  a  .5  km  diameter  cluster  at  the  TechSat  21  altitude  of  800  km.  Most  of  the  errors  from 
the  linearizations  in  Hill’s  equations  were  found  to  be  periodic  in  nature,  with  an  amplitude 
of  about  one  sixth  of  the  differential  J2©  effects.  Kong,  Miller  and  Sedwick  [22]  went  on 
to  use  these  enhanced  modelling  equations  to  analyze  minimum  propellant  configurations 
for  optical  interferometry  formations  imaging  the  Earth  from  geosynchronous  altitudes. 

Sabol,  Burns,  and  McLaughlin  [36]  used  the  enhanced  Draper  Semi-analytic  Satellite 
Theory  (DSST)  to  propagate  averaged  equations  of  motion  for  a  formation  in  the  presence 
of  a  21st  order  gravitational  field,  lunar  and  solar  third  body,  atmospheric  drag,  and  solar 
wind  perturbations.  Their  work  cataloged  the  basic  geometries  of  symmetric  formations, 
from  the  simplest  in-plane  formation  to  projected  circular  orbits,  all  with  an  eye  towards 
required  propellant  usage  for  formation  station  keeping. 

2.2  Wiesel-Floquet  Dynamics 

The  fidelity  of  modelling  formation  dynamics  was  advanced  by  more  than  three  orders 
of  magnitude  by  the  work  of  Wiesel  [54] .  Wiesel  recognized  that  the  elements  of  a  formation 
follow  periodic  paths  about  the  local  origin,  and  though  the  locations  change  over  time, 
they  are  bounded  (assuming  no  secular  growth).  He  devised  a  method  for  translating  an 
inertial  frame  Hamiltonian  to  the  same  reference  frame  used  by  Clohessy,  Wiltshire,  and 
Hill.  Keeping  the  Hamiltonian  terms  out  to  second  order  yields  a  matrix  linear  system  of 
the  familiar  form  $(f,  0)  =  A(t)$(f,  0).  Using  the  techniques  of  Floquet  [25,  29],  the  state 
transition  matrix  may  be  factored  and  solved,  see  Appendix  C. 

$(t,0)  =  F(t)eJtF(0)_1  (2.3) 

The  result  is  a  very  high  fidelity  propagation  algorithm  that  includes  potentially  all  of 
the  zonal  harmonics  and  any  other  conservative  and  non-conservative  perturbations  whose 
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potential  functions  may  be  added  into  the  Hamiltonian. 


X(t)  =  F(t)eJtF(0)~1df(0)  (2.4) 

The  thesis  by  Bordner  [1]  used  this  technique  to  generate  extremely  accurate  estimates 
of  satellite  positions  in  a  formation.  Bordner  implemented  Wiesel’s  solution  in  a  Bayes 
filter  estimation  algorithm  using  noise-corrupted  GPS  measurement  data,  carrier-phase 
differential  GPS  data,  and  inter-satellite  measurement  data. 

Wiesel  further  investigated  including  second  order  2-body  terms  and  zonal  pertur¬ 
bations  in  the  Floquet  solution  [52],  He  found  that  sectoral  and  tesseral  harmonics  only 
perturb  the  periodic  orbit,  leaving  the  relative  motion  unaffected.  The  second  order  2- 
body  and  zonal  terms  do  modify  the  F(t)  matrix,  but  do  not  introduce  instabilities  into 
the  relative  motion  dynamics.  This  formulation  also  introduced  a  change  of  variables  to  a 
set  of  six  modal  variables  Z.  see  section  C.3.  The  physical  interpretation  of  these  modal 
variables  is  more  intuitive  than  the  raw  X(t)  vector  from  the  first  formulation  [54],  and 
yields  a  dynamics  description  in  the  familiar  linear,  constant  coefficient  form 

^-Z  =  J  Z.  (2.5) 

dt  K  ’ 

The  solution  further  reveals  two  constraints,  requiring  two  of  the  modal  variables  to  be 
0  to  prevent  secular  growth  in  the  relative  satellite  paths  (yielding  the  four  degrees  of 
freedom  cited  by  Schaub  and  Alfriend).  In  a  third  advancement  [51],  Wiesel  used  his 
higher  fidelity  modelling  to  determine  a  theoretically  feasible  long-term  control  usage  level 
for  station-keeping  within  the  cluster.  The  level  of  required  control  usage  falls  easily  within 
the  capabilities  of  current  proposed  spacecraft  designs. 

2.3  Summary 

The  models  available  for  describing  relative  satellite  motion  from  the  simple  Clohessy- 
Wiltshire  solution  to  the  Wiesel-Floquet  dynamics  vary  greatly  in  sophistication  and  fi¬ 
delity.  Each,  however,  reveals  a  maximum  of  four  degrees  of  freedom  in  choosing  the 
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configuration  of  a  satellite  formation.  Adding  an  arbitrary  assumption  of  planar  forma¬ 
tions  reduces  this  to  two  since  the  assumption  imposes  two  more  constraints. 
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III.  Interferometry  Basics 

Optimizing  a  satellite  cluster  for  the  specific  mission  of  interferometry  requires  a  more 
than  cursory  knowledge  of  the  science  and  techniques  of  imaging  at  radio  wavelengths 
via  distributed  receivers.  This  chapter  begins  with  the  earliest  construction  by  Michelson 
of  a  one-dinrensional  optical  interferometer,  then  expands  the  description  to  incorporate 
multiple  receivers.  Modelling  the  aggregate  performance  of  multiple  receivers  demonstrates 
how  the  target’s  intensity  map  is  imbedded  within  a  two-dimensional  Fourier  transform  of 
the  image.  The  necessary  frequency  and  density  of  measurements  of  this  transformed  image 
(spatial  frequency  measurements)  are  described  to  yield  insight  into  nature  of  “important” 
and  “redundant”  baselines,  and  the  spatial  frequencies  they  measure. 


3.1  Michelson  Interferometry 

The  use  of  interferometric  techniques  date  back  to  the  work  of  Albert  Michelson1  who 
applied  them  to  determine  the  diameter  of  some  of  the  nearer  and  brighter  stars  (Arcturus 
and  Betelgeuse).  The  basic  interferometer  combines  two  samples  of  incident  radiation 
from  a  distant  source  into  one  receiver.  The  image  plane  of  an  optical  telescope  forms  a 
rudimentary  correlator.  In  figure  3.1,  incident  light  a\  and  b\  travel  equal  distances  to  the 
receiver.  Light  02  and  62,  approaching  at  an  incident  angle  9,  travel  paths  that  differ  by  a 
distance,  A  =  Dsin9,  and  are  therefore  out  of  phase,  \F,  by 

*=(l)Dsin  M 


The  composite  image  from  the  two  signals  displays  light  and  dark  patterns  due  to  inter¬ 
ference.  Michelson  proposed  that  the  quality  of  these  fringes  could  be  quantified  by  the 
fringe  visibility  function,  Vfr\nge,  in  terms  of  the  intensity  of  the  light  and  dark  bands,  /max 
and  /min) 


Vf 


fringe 


^max  In 
^max  +  In 


(3.2) 


1This  is  the  same  Michelson  interferometer  later  used  for  the  first  accurate  measurement  of  the  speed  of 
light  in  1879,  and  in  the  famous  Michelson-Morley  experiment  that  ended  ether-based  theory,  eventually 
leading  to  Einstein’s  relativity  theory. 
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Figure  3.1  The  basic  Michelson  Interferometer  combines  reflected  radiation  from  two 
mirrors  into  a  single  receiver.  The  constructive  and  destructive  interference 
generates  a  fringe  pattern,  a  simple  function  of  the  angle  0  and  the  baseline 
distance,  D  (or  when  measured  in  wavelengths,  D\  [46]). 

The  period  of  these  fringes  is  a  direct  measure  of  the  small  angle  0  and  the  wavelength  of 
the  radiation.  This  rudimentary  interferometer  can  measure  the  small  angle  0  by  varying 
the  interferometer  baseline,  D.  As  the  baseline  is  transformed  into  different  domains,  it  is 
nearly  always  expressed  as  a  unitless  number  of  wavelengths,  therefore  Ac  is  converted  to 
a  subscript  on  the  baseline,  or  D\= 

If  the  width  of  the  star  is  small  compared  to  0,  the  entire  image  is  covered  with  light 
and  dark  bands.  The  dark  bands  arise  as  the  phase  difference  becomes  7 r. 

7T 

0 

For  stars  comparable  in  width  to  the  spacing  between  maxima,  the  image  can  be  thought 
of  as  resulting  from  a  continuous  series  of  points  across  the  stellar  disk.  The  maxima 
and  minima  from  different  points  do  not  coincide,  and  the  images  blend  into  a  smooth 
continuum.  Two  interfered  images  of  a  circular  source  have  a  Ffringe  determined  [13]  by  a 


^  )  D0 


2D 


(3.3) 
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first  order  Bessel  function. 


DnO 


By  tuning  the  interferometer  until 
the  image  is  termed  “resolved.”2 

=  3.83.  Therefore  a  resolved 
to  satisfy 

0  =  1.22^.  (3.5) 

The  interferometer  may  be  used  to  resolve  smaller  sources  by  increasing  D  (therefore  the 
phase  difference  for  smaller  9)  or  using  shorter  wavelengths,  A0. 


Vf 


fringe 


=  2 


Dn9 
A  0 


(3.4) 


the  image  is  a  single  maxima  with  no  fringes,  Vfringe  =  0, 
The  first  order  Bessel  function  has  its  first  J\  =  0  at 
image  first  occurs  when  the  baseline  is  sufficiently  long 


3.2  Multi-receiver  Interferometers 


Modern  interferometers  include  far  more  than  two  separated  receivers,  often  on  move- 
able  supports.  The  Very  Large  Array  (VLA)  on  the  Plains  of  St.  Agustin,  west  of  Socorro, 
New  Mexico,  has  27  individual  receivers,  each  capable  of  limited  repositioning. 


Figure  3.2  The  Very  Large  Array  (left)  near  Socorro,  New  Mexico  has  its  27  antennas 
mounted  on  rails  for  limited  repositioning  capability  to  72  different  antenna 
stations.  The  Nobeyama  Solar  Radio  Observatory  (right)  in  Nobeyama, 
Japan  has  its  84  antennas  permanently  spaced  for  optimum  observations 
of  the  solar  disk. 


2Betelgeuse,  the  first  star  resolved  by  Michelson’s  interferometer  is  a  known  distance  away.  Using  the 
measured  9,  its  diameter  was  calculated  to  be  roughly  280  times  that  of  our  Sun. 
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3.2.1  Coherence,  Correlation,  and  the  van  Cittert-Zernike  Theorem.  The  simple 
Michelson  interferometer  was  adequate  for  determining  a  one-dimensional  measurement: 
diameter  of  the  source.  It  is  not,  however,  sophisticated  enough  to  yield  an  image  equiv¬ 
alent  to  that  obtained  by  an  instrument  with  a  single  large  aperture.  The  technique  for 
doing  so  requires  multiple  measurements  of  interfered  signals,  measurements  of  the  signals’ 
coherence. 

To  see  how  these  coherence  measurements  arise,  begin  by  defining  the  generic  coher¬ 
ence,  T,  of  a  monochromatic  waveform,  of  field  strength  U,  received  at  two  points  Pi  and 
P2  at  times  differing  by  r  by 

r(Pi,  P2,  t)  =  lim  [T  U(Put)U*(P2,t  +  T)dt  (3.6a) 

T— >00  Z  1  J —X 

=  (U(P1,t),U(P2,t  +  T)).  (3.6b) 


The  intensity  of  the  waveform  is  just  its  self  coherence  at  any  point,  I(P)  =  T(P,  P,  0),  or 
the  square  of  the  field  strength  [13]  J(P)  =  \U(P)\2.  Since  the  generic  wave  function  for  a 
wave  propagating  in  the  z  direction  is  U(P,t)  =  Ue^~z~ut\  the  coherence  function  at  Pi 
and  P2  may  be  determined  by  direct  calculation  as 


r(Pi,P2,r) 

r(Pi,P2,r) 


hm  —  [T 

T^oo  2 T  J_T 

lim  —  fT  e*[x(2i -*a)-«t+wt4«r]dt 

T^oo  2 T  J_T 

|f/|2e*[x(2i-22)+H. 


dt. 


Now  consider  two  incoherent  point  sources,  both  being  received  at  Pi  and  p>.  as 
in  figure  3.3.  By  the  superposition  principle,  the  received  field  strength  at  each  point  is 
simply  U  =  Ua  +  Uj,.  Since  the  discrete  point  sources  are  not  strictly  on  s,  these  functions 
no  longer  propagate  strictly  along  z.  The  component  of  propagation  may  be  calculated 
by  the  inner  product  of  the  look  vectors  from  each  receiver,  sa  and  and  their  location 
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Figure  3.3  Adapted  from  Rohlfs,  since  the  field  strengths  of  two  incoherent  sources  are 
uncorrelated.  Even  when  measured  at  the  same  point,  their  coherence  func¬ 
tion  is  0. 

vectors,  f\  and  f2-  Dropping  the  phase  offset,  their  wave  functions  become 

Ua  =  Uae^^ 

(3-8) 

Ub  =  Ube(- f 

Their  coherence  function  is  almost  as  simple. 

(U(P \,h)U*(P2,t2))  =  {[Ua(Pu  h)  +  Ub(Pu  h)},  [ Ua(P2 ,  t2)  +  Ub(P2,t2)}) 

=  h),  ua(p2,h ))  +  (ua(put  i),  ub(P2,  t2))  + 

{Ub(P1,h),Ua(P2,t2))  +  (Ub(Pi1t1),Ub(P2lt2)) 

Since  the  sources  are  incoherent,  the  field  strengths  Ua  and  Ub  are  uncorrelated,  so  the 
cross  terms  are  zero  and 


r(Pl,P2)  =  (Ua(Pl,tl),Ua(P2,t2))  +  (Ub(Phtl),Ub(P2,t2)) 

=  I  +  |Pfe|2e(^'(rW2))). 
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Extending  this  to  M  sources,  the  coherence  function  is  simply  a  function  of  the  distance 
between  the  receivers,  D  =  r\  —  r2- 

M 

T(D)  =  | C/n(s)|2e(^(s"”'D)) 

n=  1 

In  the  limit  as  M  — >  oo,  simply  integrate  across  the  entire  area  of  the  source  (the  source 
goes  to  zero  outside  a  small  angle,  so  the  limits  may  be  set  to  ±oo  without  loss  of  generality) 

/OO  /*oo 

/  |C/(s)|2e(^(s"'D))dO  (3.12) 

-OO  J  —  OO 

where,  U(s)  is  the  field  strength  per  solid  angle  (steradian).  Note  that  |17(s)|2  is  just  the 
intensity  of  the  source  at  each  point  s  =  (£,77),  so  |C/(s)|2  =  I(<^,  77) .  The  dot  product 
s  •  D  may  be  simplified  as  well.  Since  it  is  a  scalar  value  derived  from  vectors  whose 
components  are  expressed  in  different  coordinate  systems,  its  value  is  not  affected  by  the 
unitary  matrix  necessary  to  rotate  from  one  to  the  other.  Therefore  it  may  be  written 
simply  as  s  ■  D  =  |s||D|  cos(ZsD)  =  +  Ayr/. 

/oo  /*oo 

/  I(£,r 7)e(ir(A^+Aw))^dr?.  (3.13) 

-00  J  —oo 

This  is  the  monochromatic  version  of  the  van  Cittert-Zernike  Theorem.  It  shows  how 
the  coherence  function  for  a  monochromatic  wave  field  that  is  built  up  from  plane  waves 
is  related  to  the  intensity  map  [35].  This  intensity  map3,  1(^,7/),  is  the  desired  image.  It 
may  be  solved  for  if  T(D)  can  be  measured  over  a  sufficiently  diverse  set  of  D.  Measuring 
interfered  signals  gives  just  such  a  set  of  coherence  function  measurements. 

3.2.2  Synthesizing  an  Aperture.  Terrestrial  interferometry  arrays  synthesize  a 
large  aperture  telescope  using  a  large  array  of  discrete  receivers.  The  signals  from  multiple 
pairs  of  receivers  are  sent  to  a  correlator  (a  device  that  multiplies  and  time-averages  the 
incoming  signals).  The  correlator’s  response  is  not  quite  a  direct  measure  of  T(D)  since 

3Texts  on  interferometry  refer  to  the  structure  of  an  image  over  the  viewing  angles  (£,  y)  as  either 
the  intensity,  I(^,  77) ,  or  the  brightness,  B(£,  77).  By  Rohlfs  [35],  the  two  terms  are  interchangeable,  but 
“intensity”  is  the  preferred  term  here. 
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antenna  beam  patterns  are  not  perfect  and  cannot  discriminate  discrete  points  (£,  rf)  in 
the  source. 

Consider  the  delay  time,  r,  for  an  inbound  signal  as  received  by  a  pair  of  receivers, 
figure  3.4,  from  a  differential  area  of  the  source,  I(£,?y).  From  Rohlfs  [35],  the  differential 


Figure  3.4  The  delay  time,  r,  between  receivers  receiving  signals  from  s  is  a  simple 
function  of  the  baseline  dot  product  with  s. 

receiver  response,  dR,  due  to  this  delay  is 

dR(D)  =  Ajv(£,  t?)I(£,  rj)e2ini,Td^dr],  (3.14) 

where  Aat(£,  ?y)  is  the  normalized  antenna  beam  pattern  specific  to  the  type  of  receiver  and 
antenna  employed.  (If  antennas  with  infinite  apertures  were  available,  A at(£,  r/)  =  1,  and 
it  would  be  possible  to  directly  measure  T(D).)  Integrating  this  response  over  the  source 
yields 

/°0  roo  D  _ 

/  AN(£,ri)I(t,ri)e(2i™^)dtdr,  (3.15a) 

-OO  J  —  OO 

/oo  roc 

/  AN{C,d)^,ri)e^(A^+Ayn^d^dr].  (3.15b) 

-OO  J  —  OO 

In  actual  systems,  well  made  antennas  have  A n{£,v)  g°ing  t°  essentially  zero  outside  of 
a  very  narrow  beam,  so  the  integration  limits  may  be  safely  set  to  infinity.  Note  that  this 
is  not  quite  the  coherence  function  of  equation  3.13,  but  it  is  close.  It  differs  only  by  the 
factor  Ajv(£,77),  a  known  function  of  the  antenna. 

With  sufficiently  extensive  measurements  of  R(D)  over  large  variations  in  D,  a  func¬ 
tion  mapping  T(D)  may  be  synthesized  with  high  fidelity.  Further,  the  inverse  transform 
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of  equation  3.13  yields  the  high  resolution  image,  I(£,  rj),  that  a  single  large  aperture  would 
have  given. 

In  practice,  the  process  of  obtaining  an  image  involves  taking  the  Fourier  transform 
of  the  interferometer’s  coherence  measurements  at  all  available  spatial  frequencies,  deter¬ 
mined  by  the  geometry  and  spacing  of  the  receiver  elements.  As  the  correlator  records 
these  interfered  waveforms  of  the  same  wavefront  received  at  various  separations,  it  yields 
spatial  frequency  measurements.  Different  measurements  can  then  be  gathered  by  chang¬ 
ing  the  baseline  separation  between  the  receivers.  For  arrays  on  the  ground,  the  Earth’s 
rotation  provides  these  changes  by  moving  the  entire  array,  rotating  the  receiving  pairs 
through  a  predictable  locus  of  viewing  geometries. 

The  \  factor  in  equation  3.15  suggests  the  convention  of  defining  receiver  locations 
at  any  time  in  terms  of  dimensionless  numbers  of  wavelengths  rather  than  units  of  length, 
figure  3.5.  Appendix  D  summarizes  how  this  local  coordinate  system  may  also  be  in- 


/l 


Figure  3.5  Applying  the  factor  j  to  both  x  and  y  coordinates  defines  a  reference  frame 
in  terms  of  dimensionless  numbers  of  wavelengths.  Its  coordinates,  (u,  v ),  are 
defined  by  the  target  center  (£0,  r/0),  and  are  the  antenna  spacing  components 
orthogonal  to  the  viewing  direction. 


terpreted  as  components  of  spatial  frequency  [46].  For  a  source  located  at  some  angles 
(£,  rj)  from  the  array  normal,  figure  3.5  define  the  values  u  =  Dcjjs and  v  =  —  1)0  as 
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the  components  of  antenna  spacing  orthogonal  to  the  view  direction  (the  w  axis).  The 
aperture  plane,  or  ( u  v)  plane,  is  just  the  Fourier  transform  of  the  image  plane  over  its 
viewing  angle  variables  (£,77).  The  baseline  distance,  D\  from  figure  3.1,  between  receivers 
in  the  (u  v  w)  reference  system  has  components  ( u,v,w ),  all  in  dimensionless  numbers 
of  wavelengths.  Figure  3.6  shows  this  basic  geometry  of  a  two-receiver  array,  imaging  a 
distributed  source,  S. 


Figure  3.6  The  spatial  frequency  coordinate  system  used  to  describe  a  multiple  receiver 
interferometer  is  the  (w  v  w)  frame.  The  frame’s  orientation  directs  the  w 
axis  toward  the  source,  S,  and  arrays  the  receivers  about  in  the  (u  v)  plane. 
The  receiver  locations  determine  spatial  frequency  measurements. 

Most  research  on  interferometry  construction  focuses  on  receiver  motion  within  the 
(u  v)  plane  since  this  directly  determines  which  spatial  frequency  measurements  will  be 
recorded  by  the  array.  To  appreciate  the  spacing  loci  required  to  image  a  source,  consider 
the  use  of  a  two-dimensional  discrete  Fourier  transform  to  image  an  object  of  size  (A£,  A?/). 
To  reconstruct  such  an  image,  measurements  are  required  on  a  rectangular  grid  [27,  46] 
with  spacing  no  greater  than  Aw  =  ^  and  Aw  =  ^ .  Unless  receivers  have  exactly 
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measured  the  coherence  function  at  each  of  these  locations  (never  the  case),  the  values 
must  be  interpolated  from  measurements  obtained  as  the  receivers  swept  through  some 
locus  in  the  plane  (see  figure  3.7).  Therefore,  it  is  obviously  desirable  to  fill  the  (u  v)  plane 


Figure  3.7  To  image  a  source  of  size  (A£,  A 77)  at  a  bandwidth  of  AA  centered  at  A0, 
measurements  in  the  (u  v)  plane  are  required  at  spacing  no  larger  than 

Au=  TFA?  and  Av=  T757- 

as  densely  as  possible  with  measurements.  Increasing  the  size  of  the  (u  v)  plane  allows 
larger  sources  to  be  imaged,  or  smaller  sources  to  be  imaged  with  greater  resolution. 

For  optical  wavelength  instruments,  most  practical  arrays  rigorously  control  their 
receivers  to  physically  remain  within  the  (u  v)  plane,  since  removing  w  components  from 
the  coherence  measurement  requires  knowledge  of  the  signal’s  phase.  Location  tolerances 
for  radio  frequencies  are  more  relaxed  since  it  is  possible  to  accurately  measure  the  phase 
of  longer  wavelengths.  This  is  currently  not  possible  for  the  short  wavelengths  of  visible 
light. 
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The  quality  and  distribution  of  each  measurement  taken  in  the  (u  v)  plane  may  be 
more  completely  specified  by  defining  the  spatial  frequency  analog  of  the  intensity  map. 
From  Thompson  [46],  rewrite  the  response  function  of  equation  3.15  in  terms  of  variables 
(£',?/)  for  integration  about  each  set  of  viewing  angles  (£,??). 


m,v) 


cos  [2nD\  cos  6]  A(£  —  £ ,  rj 


rf)I  tf,Tf)dgdrf 


(3.16) 


This  allows  the  response  function  or  measurement  to  be  expressed  as  a  more  general  visi¬ 
bility  function  of  the  spatial  frequency  variables,  V(u,  v,  w)  by  the  general  two  dimensional 
Fourier  transform4 


V(u,v,w)  = 


Ljv(C,  r,)  I(e, 


d£drj 


'  —  OO  J  —  OO 


V7!-^2- 


(3.17) 


This  function  is  truly  a  measure  of  how  visible  each  spatial  frequency  measurement  is.  If  the 
viewed  source  (A£  and  A  rj)  is  small  enough,  then  ^  \/l  —  £2  —  ?y2  —  1^  w  ~  —  \  (£2  +  r f)  w 
and  the  visibility  function  is  simplified  in  the  (u  v)  plane  [6]  to  V(u,v)  by 


V(u,v)=  r  r  A N(tvm,v)e-2"^t+m)-r=  ■■  (3.18) 

J- OO  J- OO  V1  -  £  -  V 

The  entire  mathematical  relationship  between  V(u,v),  I (£,77),  the  received  power,  and 
the  spectral  sensitivity  are  summarized  in  figure  3.8'5.  Note  how  in  the  (u  v)  domain, 
the  spectral  sensitivity,  W(u,v),  is  just  the  transfer  function  of  linear  system.  In  control 
theory,  a  linear  system’s  transfer  function  is  the  transform  of  the  response  of  a  system  to  a 
temporal  delta  function  (which  transforms  to  a  uniform  input  across  all  frequencies).  An 
interferometer’s  transfer  function  is  the  two-dimensional  transform  of  the  system’s  response 


4Rohlfs  [35],  who  actually  cites  Thompson  frequently  in  her  work,  defines  the  transformation  from  the 
(£,  ff)  domain  to  the  (u  v)  domain  with  an  e2l,r  and  the  inverse  transformation  with  e~2ln .  Either  convention 
will  work,  provided  it  is  applied  consistently.  This  work  adopts  Thompson’s  convention. 

sThough  this  figure  is  directly  from  Thompson,  it  displays  a  minor  inconsistency.  Equation  3.18,  also 
from  Thompson,  obviously  includes  the  effects  of  the  antenna  beam  patterns  in  what  is  labelled  the  “Vis¬ 
ibility  Function.”  Equation  3.18  apparently  describes  the  Fourier  transform  of  the  “Output  Map,”  or  the 
“Measured  Visibility.” 
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Figure  3.8  From  Thompson  [46],  the  Spectral  Sensitivity,  Visibility  Function,  Power 
Reception  Pattern,  and  Brightness  Distribution  are  all  related  by  Fourier 
transformations. 


to  a  two-dimensional  delta  function  in  the  spatial  domain  (which  transforms  to  a  uniform 
input  across  all  spatial  frequencies).  A  true  point  source  would  generate  a  perfect  wave 
front,  and  thus  be  mapped  to  a  uniform  input  at  any  spatial  frequency.  The  frequency 
domain  for  the  interferometer  transfer  function  is  one  of  spatial  frequency  ( radian )  •  The 
W(u,  v )  is  a  direct  description  of  how  sensitive  the  system  actually  is  to  these  frequencies, 
hence  its  other  name,  “spectral  sensitivity.” 

The  terms  “high”  and  “low”  spatial  frequencies  may  seem  unfamiliar.  Their  contri¬ 
butions  to  the  final  image  quality  are  readily  illustrated  by  qualitatively  separating  each  in 
the  construction  of  an  example  image,  figure  3.96.  By  themselves,  the  low  spatial  frequen¬ 
cies  build  the  relatively  poor  image  resembling  that  offered  by  a  small  aperture  system. 
The  high  spatial  frequencies  refine  the  small  details  (rapid  changes  in  the  spatial  domain, 

®This  series  of  images  was  created  using  the  Filter  Design  and  Analysis  Toolbox  tool  in  MatLab’s  Signal 
Processing  Toolbox  5.0. 
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Figure  3.9  The  overall  image  quality  is  a  composite  of  the  regional  intensities  offered  by 
low  spatial  frequencies  and  the  edge-defining  properties  of  the  high  spatial 
frequencies. 

often  called  “edge  detection”)  but  are  insufficient  by  themselves  to  build  the  complete 
image. 

Historically,  the  layout  of  terrestrial  based  interferometer  arrays  sought  to  optimize 
high  and  low  spatial  frequency  measurements  and  minimize  interpolation  by  sampling 
V(u,v)  as  uniformly  as  possible,  through  all  available  spatial  frequencies.  Evenly  spacing 
the  loci  traced  out  by  the  receivers  in  the  (u  v)  plane  do  minimize  the  interpolation 
required  to  map  out  the  coherence  function  over  the  plane.  Redundant  loci  and  redundant 
spacings  add  little  to  the  resolution,  and  were  therefore  avoided  for  computational  (and 
fiscal)  frugality.  This  drove  early  arrays  to  seek  minimum  spacing  redundancy.  An  early 

^  ^  I  I  1^1  I 

Figure  3.10  This  five  antenna  array  suggested  by  Bracewell  only  has  redundancy  in  the 
unit  spacing. 

paper  by  Bracewell  [2]  presented  a  1-dimensional  spacing  scheme  for  up  to  5  receivers  in  a 
discretized  linear  array  with  minimal  redundancy  (figure  3.10).  More  capable  arrays  added 
rail  systems,  allowing  the  antenna  spacing  to  be  varied  even  more.  Using  these  arrays  to 
image  stellar  sources  entails  observation  times  on  the  order  of  12  hours.  Half  of  one  Earth 
rotation  traces  out  a  locus  for  the  receiver  pairs  that  just  repeats  during  any  portion  of 
the  next  half  rotation  (with  the  receiver  positions  reversed). 

The  “Y”  shape  of  the  VLA,  figure  3.11,  shows  some  of  the  linear  geometry  dictated 
by  the  limited  railroad  hardware,  but  is  also  the  result  of  a  0th  order  optimization  study.  As 
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Figure  3.11  The  “Y”  shape  of  the  VLA  was  chosen  to  minimize  the  number  of  (u  v ) 
cells  not  intersected  by  a  receiver  locus  over  a  viewing  period. 

described  by  Napier  [27],  Mather  discretized  the  (u  v )  plane  into  viewing  “cells”  and  com¬ 
puted  the  27  antenna  loci  over  a  viewing  period  of  ±4  hours.  The  optimization  minimized 
the  number  of  viewing  cells  not  intersected  by  one  of  the  receiver  loci  through  one  viewing 

(e?7\  071 

2 )  =  =  351 

antenna  pairs  at  any  time.  Each  antenna  pair  provides  two  measures  of  fringe  visibility, 
which  are  complex  conjugates  of  each  other,  symmetric  about  the  (u  v)  origin. 

The  702  points  actually  constitute  an  instantaneous  transfer  function,  W(u,v),  of 
the  array  [45]  (the  702  discrete  spatial  frequency  responses  to  an  impulse  function,  often 
referred  to  as  the  “snapshot”  mode  of  the  VLA).  These  are  concentrated  along  six  radial 
lines  resulting  from  the  three  linear  lines  of  antennas. 

3.3  Summary 

The  science  of  interferometry  is  concerned  principally  with  positioning  pairs  of  re¬ 
ceivers  to  measure  received  signals  over  enough  spatial  frequencies  to  yield  an  image  re¬ 
sembling  a  large  aperture’s  rendering  of  the  source.  Since  the  relative  receiver  locations 
give  measurements  that  are  Fourier  transforms  of  the  image  frequencies,  the  techniques 
naturally  address  aperture  locations  in  terms  of  the  spatial  frequencies  measured.  From 
the  earliest  interferometer  that  measured  a  single  spatial  frequency,  the  science  has  grown 
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to  large  arrays.  Increasing  the  number  of  individual  receivers  in  an 
the  number  of  measurable  baselines  at  any  instant,  since  receivers1 
pairwise.  Moving  individual  receivers  or  the  entire  array  offers  still 
me  asur  ement  s . 


array  rapidly  increases 
outputs  are  interfered 
more  spatial  frequency 
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IV.  Satellite- Based  Interferometry 

This  chapter  addresses  moving  interferometry-based  receivers  to  orbital  platforms  as  a 
logical  extension  from  terrestrial  based  systems.  Orbiting  receivers  offer  both  freedom 
from  atmospheric  distortion  and  the  potential  for  extremely  long  baselines.  This  promise 
has  prompted  much  previous  work,  and  has  resulted  in  the  basic  formulation  of  simple 
metrics.  These  metrics  as  developed  for  simple  symmetric  formations,  are  disassembled 
and  addressed  piece  by  piece,  then  generalized  to  include  time  variations. 

4.1  Previous  Proposals  and  Investigations 

Satellite-based  interferometry  is  conventionally  divided  into  the  two  regimes  of  optical 
and  radio  systems.  This  division  is  principally  driven  by  the  need  for  temporal  coherence  in 
the  received  signals.  Optical  systems  typically  preserve  this  coherence  by  tightly  controlling 
the  optical  pathways  from  receivers  to  correlator  forcing  them  to  be  identical.  For  high 
contrast  fringes,  this  equates  to  physically  controlling  pathlengths  down  to  the  order  of 
2t,  where  A0  is  the  mean  design  wavelength  and  A  A  is  the  passband  for  the  system  [42] 1 . 

The  short  wavelengths  of  visible  (and  infrared)  light  nearly  always  require  practical 
designs  to  incorporate  some  physical  connection  between  the  receivers  or  mirrors  to  rigidly 
control  the  optical  pathlengths  [28] .  Radio  wavelength  systems  have  a  more  relaxed  toler¬ 
ance  on  the  receivers.  Thus  the  receivers  may  be  physically  disconnected,  provided  that 
their  relative  positions  can  be  determined  to  less  than  the  proscribed  pathlength  tolerance. 

As  early  as  1983,  Preston  et  al.  [33],  seeing  the  capabilities  of  the  newly  completed 
VLA,  suggested  placing  at  least  one  element  of  a  terrestrial  array  into  near-Earth  orbit  and 
using  Very  Long  Baseline  Interferometry  (VLBI)  techniques  to  correlate  its  signals  with 
those  received  on  the  ground.  Such  an  orbiting  element  would  not  only  greatly  increase  the 
array’s  baseline,  but  offer  improved  (u  v)  coverage  of  low  declination  sources,  and  shorten 
required  viewing  times  (well  below  the  12-hour  dwell  times  typical  of  surface  arrays)  for 
short  duration  targets.  The  coverage  of  such  a  system  was  nowhere  near  uniform,  and  relied 

1Stachnik  and  Gezari’s  actual  tolerance  is  given  as  A^AA,  but  a  quick  dimensional  analysis  shows  this 
to  be  an  apparent  typographical  error. 
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heavily  on  the  terrestrial  elements  to  fill  in  the  center  of  the  (u  v)  plane.  One  suggested 
method  of  reducing  the  {u  v )  “holes”  was  simply  adding  more  orbiting  receivers,  essentially 
driving  planners  to  build  a  constellation  (with  minimal  reconfiguration  capability,  and  some 
very  long  wait  times,  ~20  days). 

DeCou  [9]  used  the  basic  C-W  dynamics  model  to  analyze  a  trio  of  vehicles  (two 
receivers,  one  correlator)  imaging  an  arbitrary  w  direction.  His  development  assumed  a 
constant  thrust  engine  and  yielded  a  scheme  to  gradually  fill  in  the  ( u  v )  plane  by  spiraling 
the  loci  outward  from  an  initial  starting  geometry  (an  optimistic  scheme  also  suggested 
by  Kibblewhite  [19]  back  in  1984).  DeCou’s  interferometer  description  was  for  optical 
wavelengths,  so  he  accordingly  insisted  that  the  receivers  be  controlled  to  remain  within 
the  (u  v)  plane  (the  required  fidelity  of  control  was  not  addressed).  He  also  developed  a 
methodology  for  repositioning  the  formation  to  image  a  different  desired  inertial  target. 

Johnston  and  Nock  [17],  from  the  Jet  Propulsion  Laboratory,  analyzed  the  require¬ 
ments  for  another  three  vehicle  formation,  intended  to  perform  optical  interferometry. 
Since  the  vehicles  were  not  physically  connected,  much  of  the  work  centered  on  ensuring 
the  equal  optical  pathlengths  between  the  two  telescope  vehicles  and  the  single  combining 
station  vehicle.  Using  the  basic  C-W  solution  (with  origin  at  the  combining  station  vehi¬ 
cle),  it  was  possible  to  generate  a  methodology  for  selecting  a  free- flying  formation  that 
kept  the  two  pathlengths  equal,  driven  by  the  location  of  the  target  source,  the  length  of 
the  baseline  (a  modest  10  kilometers,  maximum),  and  the  angular  velocity  of  the  reference 
orbit.  Like  DeCou,  maintaining  the  optical  pathlengths  was  all  important,  with  no  metrics 
devised  for  image  quality. 

The  notion  of  flying  an  optical  interferometer  was  formalized  in  the  Spacecraft  Ar¬ 
ray  for  Michelson  Spatial  Interferometry  (SAMSI)  program  [42] .  This  system  architecture 
consisted  of  three  or  two  vehicles  (by  perhaps  combining  one  of  the  telescopes  with  the 
combiner  station),  figure  4.1.  Though  not  physically  connected,  the  precise  receiver  lo¬ 
cation  criteria  was  proposed  to  be  handled  by  laser-locating  each  telescope  vehicle,  and 
physically  contracting/extending  the  optical  pathlength  within  the  combiner  vehicle.  The 
physical  size  of  the  combiner  therefore  proved  to  be  a  limiting  factor  on  the  allowable 
variance  in  relative  orbit  location.  Filling  in  the  ( u  v )  plane  with  adequate  samplings  of 
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Figure  4.1  The  SAMSI  program  sought  to  perform  optical  interferometry  via  two  or 
three  vehicles  with  a  modest  baseline.  After  accurately  measuring  the  receiver 
position(s),  the  optical  pathlengths  could  be  physically  matched  aboard  the 
combiner. 


V(u,  v )  also  drove  the  designers  to  include  ion  thrusters  to  drive  the  receivers  around  in  an 


expanding  ( u  v )  spiral.  Sampling  criteria  for  the  short  wavelengths  limit  such  a  spiral  to 


a  diameter  of  roughly  100  meters. 


4-2  Formation- Based  Interferometry 

Employing  a  satellite  formation  for  performing  interferometry  missions  offers  several 
advantages  over  a  constellation,  a  pair,  or  a  trio.  Since  a  cluster  can  more  densely  populate 
the  ( u  v )  plane  with  receivers,  it  offers  shorter  required  observation  times  to  yield  adequate 
samplings  of  the  visibility  function  over  the  image  plane.  These  samples  may  cross  the 
complete  image  plane  rather  than  some  locus  which  excludes  the  center  (constellations 
cannot  have  receivers  passing  through  the  center  of  the  inertial  orbit  plane).  Formation 
dynamics  also  offer  low-control  use  reconfigurations  as  mentioned  in  section  2.1.  This 
allows  relatively  cheap  reorientation  for  optimal  geometries  to  view  new  targets. 

Several  recent  works  have  analyzed  potential  uses  for  formations  capable  of  receiving 
and  correlating  visible  or  radio  signals.  Each  was  faced  with  the  difficult  task  of  defining 
a  metric  for  a  desirable  quality  in  the  interferometer. 
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Mallory  et  al.  [23]  analyzed  a  formation,  actually  a  constellation  of  geosynchronous 
orbits,  selected  to  have  overlapping  ground  traces.  Since  the  ground  traces  of  geosynch 
orbits  close  on  themselves,  the  orbital  paths  naturally  filled  in  the  (u  v)  plane  as  mapped 
onto  an  orthogonal  projection  about  the  constellation’s  nadir.  Their  choice  of  a  perfor¬ 
mance  metric  was  simply  a  statistical  accounting  of  the  number  of  receiver  pairings  in  a 
discretized  (u  v)  plane  covered  over  time,  compared  to  the  total  number  of  possible  pairings 
(similar  to  Mather’s  work  on  optimizing  the  VLA  spacings).  Though  simple,  this  metric 
did  acknowledge  time  dependency. 

Sedwick,  Kong,  and  Miller  [39]  used  the  basic  C-W  solution  parameterized  into  Hill’s 
equations,  and  added  in  perturbations  from  atmospheric  drag,  solar  wind,  and  magnetic 
forces.  Their  performance  measure  of  merit  was  based  on  the  point  spread  function,  PSF. 
This  PSF  was  refined  and  integrated  by  Sedwick,  Hacker,  and  Miller  [38].  Adopting 
Thompson’s  notation  of  £  and  r]  for  image  angular  coordinates,  it  may  be  expressed 


PSF(£,  rj) 


Gain  ( D\ ,  6) 


N 


e -(2iir(uk£+VkV)) 

k= 1 


2 


(4.1) 


where  u &  and  are  locations  in  the  ( u  v  iij)  reference  frame.  The  PSF  shows  the  effects  of 
including  or  excluding  available  baselines  from  the  compilation  of  an  image.  By  discretizing 
all  possible  positions  of  satellites  in  a  plane,  Sedwick  and  Kong  created  a  reference  array 
and  performed  a  point  by  point  comparison  of  possible  formations  to  this  fully  filled  array. 
Though  image  quality  was  mentioned,  the  work  focused  on  evaluating  synthetic  aperture 
radars  for  moving  target  identification  (MTI).  The  integrated,  or  cumulative  distribution 
of  the  Probability  of  Detection,  yields  the  Availability  of  Detection  as  a  figure  of  merit. 
The  Availability  of  Detection  offers  a  measure  of  the  fraction  of  time  that  a  given  MTI 
performance  can  be  expected.  In  this  use,  the  statistics’  dependence  on  time  was  apparent, 
but  was  not  judged  to  be  critical  due  to  the  short  duration  of  MTI  missions. 

Kong  and  Miller  [20]  had  earlier  developed  an  instantaneous  performance  metric 
of  image  quality  for  an  array  restricted  to  symmetric  (u,v)  and  (—u,—v)  receiver  pair 
locations.  They  assumed  a  “stop  and  stare”  mode  (a  necessary  simplification  for  their 
work  at  the  time),  and  applied  it  to  optical  sensing.  The  intensity  map  of  a  given  source 
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may  be  transformed  to  that  received  by  such  a  restricted  pair  by  a  technique  described  by 
Hecht  [13].  Hecht  begins  by  convolving  a  two-dimensional  point  source  at  (u,v)  with  the 
Point  Spread  Function. 

/OO  p  OO 

/  8{x'  -  u)5(y'  -  v)PSF(x',y',^T])dx'dy'  (4.2) 

-OO  J  —  OO 

Kong  and  Miller  used  as  the  PSF  the  two  variable  Fourier  transform  from  equation  4.1 
with  a  gain  specified  by  the  generic  directivity,  T>(9),  of  an  optical  system.  For  a  target 
at  off-axis  angle  9,  this  incorporates  the  fringe  visibility  function  (converting  Kong  and 
Miller’s  original  notation  with  ^  =  D\)  to  yield 

+  (4-3) 

Using  this  as  the  PSF  gain,  I(£,  rj)  becomes 

I(£,  77)  =  [  [  5{x  —  u)5(y'— v)  7r(l  +  cos  9)D\  ( '^  e ~2l7r(xkt+ykv)  dx'dy'. 

J-oc  J-oo  L  V  Dxttsui9  )\ 

(4.4) 

This  relation  may  be  used  to  map  the  symmetric  pair  assumed  by  Kong  and  Miller.  Con¬ 
volving  the  impulse  functions  with  the  receiver  locations  yields 

/oo  p  00  2 

/  Six'  -  u)8{y'  -  v)V  e-^ix'S+v'v)  +  e -2in(-x'Z-y'v)  dx'dy\  (4.5) 

-OO  J  —  OO 

By  the  sifting  property  of  the  5  function  this  is  easily  seen  to  be2 

2 

I(£ ,  rj)  =  V  e~2iwiui+vv)  _)_  e-2in(-u£,-vri) 

=  V  |cos  2tt{u £  +  vrj)  —  i  sin  2tt{u^  +  vrj)  +  cos  2ir{u£>  +  vrj)  +  isin27r(ti£  +  vrj)  |2 

=  7r(l  +  cos  9)D\  (  —  f^A7r  ^  2  cos  2n{u^  +  vy) 

\  D\tt  sin  9  J 

2Kong  and  Miller’s  actual  presented  expression  was  |V(1  +  cos  6)D\  ^  2 cos7r(ig;  +  , 

but  that  result  is  respectfully  disputed  here. 
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Using  the  property  of  the  Bessel  function  that  lima_>o  Jla ^  =  \  simplifies  I (£,77).  For 
0  =  0,  the  condition  that  occurs  when  w  is  defined  to  point  directly  at  the  source, 


7r(l  +  1)D\  2  (cos  +  vrjj) 

[27rl?^  cos  2tt(u^  +  vq)]2  . 


(4.7a) 


The  contributions  of  N  symmetrically  spaced  receiver  pairs  can  be  summed  up  to  yield  an 
effective  intensity  map,  I(^,q),  via 


(27tT>a)2 

N 


'  N 

Y  (cos  27r(iifc£  +  vkq))2 
Jk= 1 


(4.8) 


which  provides  an  instantaneous  assessment  of  the  received  intensity  for  any  location  (£,  q) 
by  all  receivers  in  the  formation.  This  relation  was  further  used  to  determine  an  I0(£,  77) 
function  that  would  be  generated  if  receivers  occupied  all  symmetric  locations  ( u ,  v)  in  a 
discretized  formation  plane.  This  instantaneous  measure  of  the  best  image  obtainable  by 
a  discretized  formation  was  analyzed  over  a  discretized  m  x  m- pixel  image  field  to  yield 
the  aggregate  mean  square  error. 


MSE 


gfaigfci  (He.,  (;,.ij))2 

m  x  m 


(4.9) 


Kong  and  Miller  discretized  the  imaging  field  into  such  an  m  x  m  grid  and  performed  a 
0th  order  optimization  through  a  simulated  annealing  technique  [16] .  Each  iteration  of  this 
technique  compared  combinations  of  possible  receiver  pairings  to  the  summation,  I0(6  ,  rlj )  > 
of  all  possible  (2601  permutations  in  their  model)  symmetric  receiver  pairings  of  (u,  v)  and 
(~u,  -v). 


4-3  Analyzing  a  Comparison  Metric 

The  logic  put  forward  by  Kong  and  Miller  of  comparing  the  intensity  map  from 
a  discrete  receiver  array  to  a  filled  aperture  of  equal  size  is  almost  inarguable.  A  filled 
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aperture  is  essentially  an  infinite  number  of  receivers.  Comparing  the  discrete  receivers’ 
performance  to  this  infinite  number  basically  picks  out  which  receivers’  outputs  are  crucial, 
and  which  are  redundant.  While  advancing  the  logic  used  in  Kong,  Miller,  and  Sedwick’s 
investigation,  there  are  several  issues  that  bear  scrutiny.  The  selection  of  an  appropriate 
I0(£,77)  has  several  possibilities.  The  calculation  of  1(^,7/)  is  somewhat  complex  for  asym¬ 
metrical  formations.  Furthermore,  the  size  of  the  image  plane  over  which  the  comparison  is 
performed  has  some  effect  on  the  calculated  optimal  configuration  and  a  huge  effect  on  the 
calculation  budget  necessary  to  determine  that  optimum.  The  metric,  its  aggregate  aper¬ 
ture,  and  the  image  plane  size  are  formalized  here.  The  reference  apertures  are  addressed 
in  chapter  V. 

Performing  a  comparison  of  a  pair  of  two-dimensional  functions  over  an  area  is  nec¬ 
essarily  a  double  integration.  The  aggregate  aperture  I(£,  77,  t)  is  a  function  of  the  receiver 
locations  at  any  time,  and  l0(£,rj)  is  determined  by  the  chosen  aperture  geometry.  Em¬ 
ploying  a  second  order  metric,  this  becomes 

rv  rtf  r_  n2 

J image  (path,  t)  =  /  /  [l(£,  ??,  t)  -  I0(£,  7?)]  d^dr].  (4.10) 

J  Vo  J  £o 

Kong,  Miller,  and  Sedwick’s  report  [21]  to  the  TechSat  21  office  included  an  outstand¬ 
ing  discussion  of  optimization  techniques  based  on  a  metric  of  comparison  with  a  filled 
aperture.  Kong’s  thesis  (a  subsection  of  their  report)  developed  a  filled  circular  aperture, 
discretized  the  image  plane  into  rn  x  m  cells,  and  performed  a  point  by  point  comparison, 
using  the  metric  put  forth  in  his  paper  with  Miller  [20],  equation  4.9.  Minimizing  such 
a  metric  provides  an  image  with  a  minimum  instantaneous  error.  Their  pioneering  work 
showed  the  advantages  of  departing  from  earlier,  geometrically  symmetric  configurations, 
to  yield  minimal  errors  over  appropriate  portions  of  the  image  plane.  They  further  inves¬ 
tigated  propellant  budgets  needed  to  reconfigure  formations  to  these  optimal  geometries. 

4-3.1  Comparing  Apertures.  Using  the  image  plane  to  compare  outputs  of  differ¬ 
ent  apertures  is  essentially  a  differencing  of  two  energy  fluxes.  Each  aperture  generates  an 
intensity  map  carrying  units  of  st^adfan  >  and  the  comparison  process  requires  calculating 
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the  flux  generated  by  each  aperture  at  locations  (£,  rj)  across  the  image  plane  and  perform¬ 
ing  the  subtraction  I(£,  rj,  t)  —  I„(£,  77).  A  well  matched  pair  of  apertures  yields  only  slight 
differences  at  each  location,  and  the  normalized  sum  in  equation  4.9  is  just  an  average  of 
small  numbers. 

To  illustrate  the  nature  of  this  comparison  sum,  consider  a  one-dinrensional  filled 
aperture,  I0(£),  compared  to  a  discrete  aperture,  I(£)  (developed  in  the  next  section)  of 
just  two  receivers.  The  filled  aperture’s  map  is  readily  generated  using  equation  5.1,  and 
equation  4.11  gives  the  aggregate  aperture’s  map.  The  driving  question  arising  from  the 

1-D  Filled  vs  2  Receivers 


Figure  4.2  Placing  two  receivers  at  the  maximum  limits  of  a  filled  1-dimensional  aperture 
yields  a  poor  approximation. 

comparison  is  where  to  place  the  two  receivers  to  most  closely  match  the  two  intensity  maps 
across  the  image  line.  In  figure  4.2  this  equates  to  minimizing  the  area  between  the  graphs 
for  the  filled  and  the  two-receiver  aperture.  From  the  Michelson  interferometer  of  figure  3.1, 
the  initial  temptation  for  an  aperture  of  width  D\  is  to  place  the  two  receivers  as  far  apart 
as  possible,  or  at  the  ends  of  the  slit  D\  apart.  This  actually  yields  a  poor  approximation 
of  the  complete  filled  aperture  since  its  measurement  of  only  the  highest  spatial  frequency 
comes  at  the  expense  of  generating  large  sidelobes  in  the  reception  pattern.  These  sidelobes 
represent  reception  angles  receiving  large  amounts  of  interference.  Moving  the  receivers 
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closer  together  reduces  the  size  of  these  sidelobes,  at  the  expense  of  widening  the  central 
peak.  Across  a  finite  image  line,  there  is  certainly  an  optimum  separation  of  the  two 
receivers  to  minimize  the  sum.  Adding  more  discrete  receivers  to  the  aperture  can  shape 
the  aggregate  intensity  map  to  look  more  and  more  like  the  filled  aperture. 

For  an  actual  two-dimensional  image  plane,  the  comparison  process  is  still  the  same, 
but  another  consideration  arises.  Whereas  a  one  dimensional  reference  aperture  has  length 
as  its  only  parameter,  a  two  dimensional  aperture  has  a  size  and  a  shape.  For  space- 
based  interferometry,  the  shape  of  the  reference  should  be  influenced  by  the  achievable 
loci  of  receiver  paths.  From  the  dynamics  descriptions  of  section  2.1  and  inspection  of 
equations  2.2,  the  relative  loci  are  elliptical,  tracing  out  ellipses  with  a  2:1  ratio  of  major  to 
minor  axes.  Ignoring  these  loci  and  choosing  a  circular  shape  for  the  comparison  aperture 
is  understandable  since  it  has  a  known  solution  (section  5.2.2)  and  matches  the  shape 
of  most  primary  mirrors,  lenses,  or  filled  antennas.  Though  mathematically  simple  this 
comparison  is  best  limited  to  either  instantaneous  use  or  to  receivers  tracing  out  circular 
loci  ( not  the  case  for  a  formation).  Trying  to  force  a  formation’s  performance  to  mimic  a 
circular  aperture  spawns  the  question  of  what  to  choose  for  the  size  of  the  circle;  choosing 
either  the  major  or  minor  axis  causes  difficulties,  figure  4.3.  A  minor  axis-sized  aperture 


Figure  4.3  Comparing  receivers  tracing  out  elliptical  loci  to  a  circular  aperture  is  dubi¬ 
ous,  whether  the  circle  is  defined  as  circumscribed  or  inscribed. 
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ignores  fully  half  of  the  higher  spatial  frequencies  measurable  by  receivers  orbiting  out  to 
the  extreme  portions  of  the  ellipse.  An  aperture  sized  to  the  major  axis  compares  the 
receivers  aggregate  intensity  map  to  an  aperture  with  areas  not  reachable  by  the  dynamics 
of  the  formation.  Therefore,  matching  the  comparison  aperture  shape  to  the  receiver  paths 
ensures  that  all  achievable  baselines  are  considered  and  that  no  potentially  useful  paths  are 
ignored.  With  the  reference  aperture  properly  shaped,  it  offers  the  optimization  routine  a 
complete  set  of  possible  baselines  and  no  impossible  ones.  The  optimization  process  then 
has  the  task  of  choosing  the  most  beneficial  baselines  to  be  sampled  from  this  set. 

4-3.2  Aggregate  Apertures.  Replacing  a  continuous  aperture  with  a  set  of  discrete 
receivers  necessarily  degrades  its  performance,  but  calculating  the  aggregate  performance 
is  not  computationally  challenging.  The  Superposition  Principle  [13]  allows  the  signals 
from  each  receiver  pair  to  be  simply  added  algebraically.  Consider  the  one-dimensional 
aggregate  aperture  in  the  last  section  with  pair  spacings  Af..  Since  each  pair  produces  the 
map  of  a  point  source,  the  aperture’s  map  is  simply  a  finite  summation  of  the  real  part3 
of  the  PSF  for  each  pair.  The  normalized  aggregate  intensity  map  from  equation  4.1  for 
N  receiver  pairs  becomes 
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^e-2inAk(, 


(4.11) 


Inspection  of  equation  4.11  reveals  that  each  receiver  pair  contributes  a  (cos)2  function 
whose  frequency  is  determined  by  the  size  of  the  spacing  between  the  pair. 

In  two  dimensions,  each  receiver  pair  contributes  a  baseline  that  supplies  the  (cos)2 
distribution  in  one  direction  on  the  image  plane,  and  a  constant  along  the  perpendicular 
direction,  figure  4.4.  Mimicking  the  filled  two-dimensional  aperture  is  therefore  a  search  for 
optimal  baseline  lengths  and  orientations  such  that  the  aggregate  intensity  map  resembles 
that  generated  by  the  filled  aperture  through  a  comparison  at  each  (£,  rj)  location. 


3Due  to  the  symmetry  of  the  PSF,  its  imaginary  part  could  just  as  well  be  chosen  and  squared  to  yield 
the  intensity  without  loss  of  generality  [13]. 
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Single  Baseline  Reception 


xi  (rad) 


Figure  4.4 


Each  baseline  contributes  a  two-dimensional  (cos)2  function  parallel  to  the 
baseline. 


As  discussed  in  section  4.2,  early  work  by  Kong  and  Miller  [20]  simplified  the  ge¬ 
ometry  of  their  search  by  assuming  formations  of  matched  pairs  of  receivers.  This  leads 
to  simplified  analytical  expressions  for  the  aggregate  intensity  map  generated  by  N  pairs 
of  receivers,  but  this  simplification  was  at  the  expense  of  inherently  creating  redundant 
baselines,  figure  4.5. 

By  the  Binomial  Theorem,  for  M  =  4  receivers,  there  are  actually  P  =  =  6 

baselines.  By  forcing  M  receivers  into  symmetric  pairs,  baselines  (4)  and  (6)  are  redundant. 
This  redundancy  is  easily  removed  by  no  longer  insisting  on  symmetric  pairs.  (In  his 
thesis  [21],  Kong’s  optimizations  allowed  full  two  degree  of  freedom  for  receiver  placement.) 


4-11 


Figure  4.5  Summing  the  contributions  from  N  =  2  pairs  overlooks  four  other  baseline 
pairings  available.  Individually  locating  M  =  4  receivers  potentially  creates 
6  unique  baselines,  therefore  6  spatial  frequency  measurements. 


This  generation  of  permutations  may  be  easily  formalized.  Let  P  be  the  number  of 
pair  permutations  that  exist  in  the  set  of  M  receivers,  and  let  be  the  mapping  that 
generates  B,  the  set  of  P  pairs  without  order.  The  right  diagram  in  figure  4.5  has  a  B  of 
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Using  this  set,  it  is  possible  to  generate  the  more  complete  permutational  aggregate  image 
intensity,  Ip(^,r/,i). 

p 

=  (4.13) 

k= 1 


Recall  from  Rohlfs  [35] ,  the  absolute  position  of  each  receiver  within  the  (u  v)  plane  is  not 
important.  Measurements  of  the  coherence  function  are  strictly  a  function  of  the  difference 
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between  the  positions4  (equation  3.13).  Therefore,  for  ease  of  visualization,  translate  the 
six  permutations  to  give  origin-centered  loci. 


Figure  4.6  The  six  permutations  available  from  four  receivers  can  be  translated  to  the 
origin,  without  changing  the  spatial  frequency  measurement. 


Therefore  the  aggregate  aperture  intensity  map  of  equation  4.11  expanded  to  two 
dimensions  must  be  summed  from  1  to  P  to  include  all  permutations.  For  the  full  planar 
map,  the  normalized  intensity  becomes 
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(4.14) 


This  is  the  normalized  instantaneous  aggregate  intensity  map. 

The  size  of  the  image  plane  over  which  the  comparison  is  performed  has  a  direct 
impact  on  the  required  computation  budget.  Doubling  the  length  of  the  axes  quadruples 
the  number  of  image  intensity  points  that  must  be  calculated  and  compared  (unless  the 
discretization  is  made  more  coarse).  Note  from  figure  4.2  that  the  Bessel  function  which 
nearly  always  occurs  in  the  I0 (£>??)  for  filled  apertures  falls  off  rapidly  outside  the  central 
peak.  Indeed,  it  is  almost  negligible  after  the  third  zero  crossing.  Looking  at  a  generic 

4Rohlfs’  equations  6.10  and  6.17. 
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function  of  the  form 


JM) 

A 


2 


(4.15) 


Response  = 


reveals  zeroes  at  A  =  [3.8317  7.0156  10.1735  ...],  figure  4.7.  Limiting  the  plane  of 


Bessel  Zeros 


Figure  4.7  Graphing  the  square  of  a  generic  Bessel  function  from  equation  4.15  out  past 
10  reveals  zeros  at  A  =  [3.8317  7.0156  10.1735]. 

comparison  to  lie  within  the  first  zero  is  overly  restrictive,  and  essentially  assumes  a  very 
sparsely  filled  image  (objects  in  the  sidelobe  viewing  locations  are  few) .  Moving  out  to  the 
second  Bessel  zero  creates  a  fairly  large  image  plane.  Moving  out  to  the  third  Bessel  zero 
creates  a  very  large  image  plane  with  most  of  the  reference  function  nearly  zero.  Such  a 
comparison  plane  assumes  a  very  densely  populated  image,  and  will  drive  an  optimized 
array  towards  small  spatial  frequencies  to  yield  these  regions  of  near  zero.  This  suggests 
performing  the  comparison  over  an  image  plane  whose  size  is  determined  by  the  aperture 
D\,  such  that 

7T  U  Dx  =  ±7.0156  .  (4.16) 


This  encompasses  the  first  natural  sidelobe  or  ring  in  the  Airy  disk,  but  does  not  devote 
expensive  computation  cycles  to  vast  areas  of  near  zero  gain.  (The  effects  of  the  regions 
between  the  second  and  third  zeros  are  not  totally  ignored.  The  corners  of  a  square  image 
plane  extend  out  into  the  small  non-zero  region  between  7.0156  and  10.1735.) 
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4-4  Mean  Formation  Metrics 


Nearly  all  interferometers  increase  the  number  of  baselines  sampled  by  dwelling  on  a 
target  for  some  finite  duration,  during  which  time  the  receivers  trace  out  paths  in  the  {u  v) 
plane.  The  aggregate  intensity  map  of  the  previous  section  only  yields  an  instantaneous 
measure  of  image  errors.  For  receivers  in  constant  independent  motion,  this  metric  may  be 
generalized  to  include  these  instantaneous  errors  averaged  over  a  finite  time.  For  a  satellite 
formation,  such  a  metric  offers  optimal  receiver  paths  such  that  imaging  over  time  from 
any  start  time  will  suffer  from  minimal  errors  over  some  finite  dwell  time.  Aside  from 
reinforcing  the  need  to  consider  path  shapes  as  described  in  section  4.3.1,  this  adds  yet 
another  integral  to  the  metric. 

Jimage(path)  =  f  I  f  [lP(£,  rj,  t)  -  I0(£,  r?)] 2  d^dr/dt  (4.17) 

J  to  J  F/o 

Since  baseline  geometries  repeat  every  half  rotation,  a  logical  dwell  time  for  a  satellite 
formation  might  well  be  one  half  orbital  period,  Discretized  over  an  m  x  rn  image  plane 
and  n  time  steps,  this  becomes 

1  n  m  m 

J  image  (path)  =  „  \  ^  ^  ^  (ip  (£i>  Vj  >  *fc)  —  Io(6,?7j))  •  (4-18) 

4  >k=lj=n=l 

Note  that  as  developed,  J image  (path)  is  the  average  of  the  instantaneous  metric  over 
a  imaging  dwell  time,  not  the  effects  of  all  baselines  over  time  as  compared  to  a  perfect 
aperture.  This  is  an  important  difference,  as  it  generates  optimal  paths  for  yielding  minimal 
image  errors  for  any  start  time. 

4-5  Summary 

Moving  interferometric  receivers  into  orbit  is  a  logical  extension  from  terrestrial  based 
arrays.  The  idea  dates  back  to  the  early  1980s  and  fostered  preliminary  work  on  how  to 
compare  the  performance  of  prospective  arrays.  The  current  best  measure  of  merit  for 
an  array  is  a  summation  of  its  image  plane  intensity  map  errors  as  compared  to  a  filled 
aperture.  This  is  a  difficult  comparison.  Assuming  symmetric  receiver  placement  yields 
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simple  analytic  expressions  for  discrete  array  performance,  but  at  the  cost  of  creating 
redundant  baselines.  Generalizing  the  mathematics  to  describe  asymmetric  arrays  and 
time-dependence  is  not  difficult,  but  keeping  track  of  the  rapidly  expanding  number  of 
baselines  with  larger  numbers  of  receivers  is  challenging.  The  aperture  used  as  a  reference 
must  also  be  selected  with  the  nature  of  receiver  motion  in  mind  if  the  metric  is  expanded 
to  include  time  variations.  A  logical  image  plane  size  for  the  comparison  may  readily  be 
found  by  observing  the  nature  of  a  filled  aperture’s  I0(£,  rj).  Several  of  these  filled  aperture 
maps  are  developed  in  the  next  chapter. 
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V.  Apertures,  Apertures,  Apertures 

A  rigorous  analysis  of  an  optimal  distributed  aperture  is  not  possible  without  first  analyzing 
the  basic  physics  of  the  finite,  filled  apertures  used  as  references.  This  chapter  reviews  the 
intensity  map  generated  by  several  regularly  shaped  apertures  receiving  electromagnetic 
radiation,  then  moves  on  to  develop  the  output  of  an  asymmetric  aperture.  Each  is  derived 
by  integrating  the  PSF  over  the  aperture  geometry. 

5.1  One  Dimensional  Apertures 

The  nature  of  a  comparison-based  image  quality  metric  is  most  easily  seen  by  first 
addressing  a  1-dimensional  continuous  aperture  (a  slit).  Though  not  particularly  practical 
for  generating  an  image1  (an  ideally  focused  image  would  only  be  ideal  in  one  dimension), 
its  performance  is  easily  displayed.  The  1-dimensional  point  spread  function  is  readily 
integrated  across  a  slit  of  width  D\. 


Figure  5.1  Integration  of  the  point  spread  function  across  a  1-dimensional  aperture, 
yields  the  intensity  map  over  the  single  angular  variable,  £. 


1The  original  Michelson  interferometer  was  indeed  a  1-dimensional  system,  but  he  sought  only  to  measure 
the  widths  of  nearby  stars  not  to  generate  images  of  them.  Therefore  he  was  only  interested  in  measuring 
the  highest  spatial  frequency  available  to  yield  the  best  possible  edge-detection  of  the  stars’  disks. 
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The  shape  of  the  intensity  mapping  is  determined  by  Da  over  all  viewing  angles  £,  and 
increases  overall  with  increased  slit  width  (easily  normalized). 


Normalized  1-D  Aperture  Response 


xi  (rad) 


Figure  5.2  The  single  dimension,  Da,  of  a  one-dinrensional  aperture  determines  the 
shape  of  the  intensity  map  over  all  look  angles  (Da  =  100  wavelengths  and 
Da  =  200  wavelengths  shown). 
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5.2  Two  Dimensional  Apertures 

Filled  two-dimensional  apertures  have  some  geometry  in  the  (u  v)  plane  covered  by 
a  continuous  receiver  of  some  sort  (the  objective  lens  or  primary  mirror  of  a  telescope 
are  common  examples).  The  intensity  map  created  by  such  an  aperture  is  generated  by 
integrating  the  two-dimensional  PSF  over  the  whole  aperture  for  each  pair  of  viewing 
angles.  For  apertures  with  simple  geometries  (circles,  squares),  analytic  solutions  exist  for 
the  image  plane  intensity  maps.  An  analytic  solution  for  an  elliptical  aperture  also  has 
been  developed2,  and  the  derivation  is  shown  in  section  5.2.3. 

5.2.1  Square  Apertures.  The  mathematically  simplest  aperture  to  integrate  is 
certainly  a  square,  since  the  point  spread  function  may  be  split  into  two  independent 
functions  along  the  £  and  i]  axes. 
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2  Based  on  extensive  literature  search,  the  author  believes  this  to  be  an  original  derivation  for  this  type 
of  aperture. 
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Mapped  over  a  discretized  image  plane,  this  generates  a  planar-symmetric  intensity  map, 
figure  5.3. 


Normalized  Square  Aperture  Response 


Figure  5.3  A  square  aperture’s  intensity  map  is  the  product  of  two  1-dinrensional  maps. 

This  map  is  generated  by  a  square  aperture  D\=  100  each  side.  The  intensity 
or  response  given  over  the  image  plane  is  depicted  as  the  height  of  the  surface. 


Though  seldom  used,  this  is  the  best  response  a  square  lens  or  primary  mirror  of  dimension 
D\  could  yield  for  a  point  source. 


5.2.2  Circular  Aperture.  As  discussed  in  section  4.3.1,  a  circular  aperture’s 
intensity  map  is  not  quite  adequate  for  use  as  a  reference  for  formation  performance.  Its 
derivation  is  worthwhile,  however,  since  it  may  be  extrapolated  to  an  elliptical  aperture.  A 
normalized  intensity  map  of  an  aperture  may  be  stepwise  built  by  determining  the  Fourier 
transform  of  the  point  spread  function  from  the  spatial  frequency  domain,  normalizing  it 
by  the  maximum  power  reception  pattern,  then  squaring  it  to  yield  an  intensity  [13]. 
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With  this  approach,  the  square  aperture  derivation  may  be  modified  to  a  more  prac¬ 
tical  circular  aperture  by  a  natural  conversion  to  polar  coordinates  in  both  the  spatial  and 
spatial  frequency  domains.  Following  a  derivation  by  Rohlfs  [35] ,  convert  aperture  positions 
(u,  v)  to  the  polar  coordinates  (r,9)  and  viewing  angles  (£,  r/)  to  the  coordinates  (p,  <J>). 


u  =  r  cos  6  £  =  pcoscj) 

v  =  r  sin  9  rj  =  p  sin  </> 


Figure  5.4  In  the  spatial  frequency  domain,  convert  aperture  positions  ( u ,  v )  to  the  polar 
coordinates  (?’,  9).  In  the  spatial  domain,  shift  viewing  directions  (£,  rj)  to  the 
coordinates  (p,  </>)• 

The  non-normalized  power  reception  pattern  is  converted  by 
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where  the  reception  pattern  dependence  on  (f>  is  immediately  removed  by  axial  symmetry. 
One  integral  may  be  removed  by  rearranging, 
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and  noting  the  definition  of  the  Bessel  function  (JQ)  of  the  first  kind  of  order  0  for  any 
value  a 


1 


Mat)  =  f-  /  eiacosedO. 

2lT  Jo 


Setting  a  =  2irpr  allows  equation  5.5  to  be  rewritten  as 


(5.6) 


f{p)  =  /  MZnpr)  r  dr. 
Jo 


(5.7) 


The  normalized  intensity  map  may  now  be  constructed  by  noting  the  maximum  value  of 
Ji  is  1.  Performing  the  normalization  and  squaring  yields 


Up)  = 

A  change  of  variables  z  =  2irpr  allows3 


/02  J0(2irpr)rdr 
Rx 

Jo2  rdr 


(5.8) 


Up)  = 


[•D\irp  j  /  \  z  dz 
JO  2-kP  2ttP 

fo  2  rdr 

A^pZ  f0DXnP  Jo(z)zdz 


2 


[DU p*  J0 


8 

* D\irp  -|  2 

JQ(z)zdz 


(5.9a) 

(5.9b) 

(5.9c) 


3This  point  in  the  derivation  corrects  a  misprint  by  Rohlfs  who  asserted  D^np  5.9c 
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Again  noting  the  property  of  Bessel  functions  [13]  that 


/  aJ0{a)da  =  /3Ji(/3),  (5.10) 

Jo 

allows  the  last  integral  to  be  replaced  by  a  Bessel  function  (Ji)  of  the  first  kind  of  order  1. 


Up) 


2 

D2  2^Dx^p)Ji(DXTrp) 

2Ji(HA7rp)'  2 

D\irp 


(5.11a) 

(5.11b) 


This  intensity  map  yields  the  classic  Airy  disk  pattern  with  a  central  peak  of  width  inversely 
proportional  to  the  aperture  diameter  D\,  figure  5.5.  Note  that  the  pattern  is  symmetric 


Normalized  Circular  Aperture  Response 


Figure  5.5  The  Airy  disk  pattern  resulting  from  a  circular  aperture  of  diameter  D\  =  100 
is  shown  here  over  a  discretized  set  of  viewing  angles. 

about  the  normal  vector  of  the  image  plane  (6  does  not  appear  in  equation  5.11b). 
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5.2.3  Elliptical  Apertures.  The  intensity  map  generated  by  an  elliptical  aper¬ 
ture  may  be  extrapolated  from  a  circular  aperture  to  yield  the  reference  described  in 
section  4.3.1  shaped  to  match  the  elliptical  loci  of  formation  elements.  Begin  with  the 
non-normalized  intensity  map  of  equation  5.3,  and  assume  an  elliptical  aperture  with  its 
semi-minor  axis,  b,  and  semi-major  axis,  a,  aligned  with  the  (u  v)  axes  respectively.  The 
conversion  to  polar  coordinates  is  essentially  the  same,  figure  5.6.  With  an  asymmetric 


Figure  5.6  An  elliptical  aperture’s  intensity  can  use  the  same  conversions  to  polar  coor¬ 
dinates  in  both  domains  as  used  for  a  circular  aperture. 

aperture,  however,  the  spatial  domain  contributions  to  the  intensity  map  are  no  longer 
symmetric  about  the  aperture  normal.  This  alters  the  limits  of  the  inner  integration  to  a 
function  of  the  angle,  9,  though  the  Jacobian  for  the  conversion  is  still  just  r. 

/(£,»?)  =  [ e2lMu^+vv)du  dv  (5.12a) 

A 

1  r2n  rR(0) 

/ellipse ( ft  4>)  =  ^Jo  Jo  e2i^cos^rcose\  dr  de  (5.12b) 

The  outer  radius  of  the  aperture,  R(9),  may  be  derived  from  the  parametric  representa¬ 
tion  [44]  of  an  ellipse 

u(9)  =  b  cos  9 
v(9)  =  a  sin  9 
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The  equation  of  the  outer  radius  is  found  readily. 


R2(u,v)  =  u2(9)  +  v2(9) 

R2(9)  =  ( b  cos 9)2  +  (a  sin#)2 
R(9)  =  y/(b  cos #)2  +  (a  sin#)2 

Inserting  this  into  equation  5.12b  gives  the  basic  polar  representation  of  the  intensity  map. 

1  t-2-K  r^/Jb  cos  8)2+(a  sin  6)2 

/ellipse (#j  <j>)=7T  /  ^(pcos *)(r cos 6)  r  dr  dQ  (5.15) 

^  Jo  Jo 

Since  the  limits  of  the  inner  radius  integral  are  not  a  constant  radius,  but  a  function  of 
the  outer  variable  of  integration,  #,  the  order  of  integration  may  not  be  interchanged  to 
build  a  Bessel  function  of  the  first  kind  (as  in  the  circular  aperture).  However,  is 

only  needed  to  generate  an  intensity  map  over  an  arbitrary  image  plane.  The  image  plane 
may  just  as  well  be  square  and  mapped  in  cartesian  coordinates,  /(£,?/).  Such  an  /(£,  ?y) 
may  be  derived  two  different  ways. 

The  most  direct  and  least  efficient  derivation  of  /(£,  rj)  is  to  numerically  integrate 
equation  5.15  over  a  discretized  aperture.  A  simple  cartesian  integrator  may  be  employed 
across  a  circumscribed  square  in  the  (u  v)  plane  if  a  logic  statement  is  added  which  drives 
the  integrand  to  zero  outside  the  elliptical  radius. 

”>  =  i  ['  £'  ^+vv>  [(£  +  i) 5 :]  *•  ^  (5'16) 

This  method  is  computationally  taxing,  since  each  point  in  the  discretized  image  plane  must 
be  calculated  via  numerical  integration  across  the  entire  aperture  ellipse4.  Nevertheless,  it 
is  certainly  feasible  and  yields  an  intensity  map  with  a  familiar  appearance,  figure  5.7.  Not 
surprisingly,  the  contraction  of  the  elliptical  aperture  along  the  u  axis  causes  an  elongation 
of  the  intensity  map  peak  and  sidelobes  along  the  £  axis. 

4Figure  5.7  required  25.5  hours  of  computation  on  a  1GHz  PC,  using  MatLab’s  standard  variable-stepsize 
double  integration  function  dblquad( .  .  . ). 
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Normalized  Elliptical  Aperture  (numerical) 


Figure  5.7  An  elliptical  aperture  may  be  laboriously  integrated  numerically  for  each  set 
of  look  angles  to  yield  an  intensity  map  I0(£,  ^ellipse- 

Equation  5.15  may  be  solved  in  a  more  elegant  fashion  by  readdressing  the  problem 
of  directly  determining  /(£,  rj)  through  an  approach  resembling  one  used  by  Johannes 
Kepler.  Kepler  solved  the  difficult  problem  of  elliptic  satellite  motion  by  mapping  the 
known  solution  for  circular  motion  about  a  circumscribed  circle  onto  the  ellipse  [49] . 

By  circumscribing  a  circle  of  diameter  D\  about  the  elliptical  aperture,  it  is  apparent 
that  the  ellipse  is  simply  a  circular  aperture  contracted  by  the  ratio  of  ^  in  the  u  direction, 
figure  5.8.  This  suggests  defining  an  eccentric  radius,  r' .  The  coordinate  transformation 
in  the  spatial  frequency  domain  is  modified  to 

v!  =  cos  6 

v'  =  v  =  r  sin  9 
or 

^u'  =  rcos9 
v'  =  r  sin# 
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Figure  5.8 


Circumscribing  a  circle  about  the  elliptical  aperture  highlights  mapping  the 

circular  aperture  response  at  r  to  the  eccentric  radius  r'  =  \J u)"  +  v2. 
The  major  axis  of  the  ellipse  is  D\. 


This  is  just  the  original  polar  coordinate  system  used  to  solve  the  circular  aperture,  mod¬ 
ified  by  the  eccentricity  of  the  ellipse.  Since  the  coordinate  transformation  from  (u,  v )  to 
(r,  9)  is  the  same  as  from  ( |f/,  v ')  to  (r,  9),  the  polar  point  spread  function  may  be  written 


PSFellipse(r,  9,  p')  =  e2^'(rcosW)  (5.18) 

in  terms  of  an  image  plane  eccentric  radius,  p',  where 

P2  =  Qe)  +(V)2  (5.19a) 

P  =  (5-19b) 

The  PSF  may  be  transformed  back  to  its  Cartesian  form,  carrying  along  the  |  contraction 
along  the  u  direction. 


P  S P ellipse  ^  ?  £5  V) 

_  e2iir(u'^+v'ri) 

(5.20a) 

PSF ellipse(^)  ^5  £5  V) 

—  e2in(,%u£+vr)) 

(5.20b) 

PSF ellipse(^)  ^5  £5  V) 

—  e2i7r(“(|?)+OT?)_ 

(5.20c) 
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The  last  bit  of  algebra  highlights  the  similarity  to  the  circular  aperture,  with  the  eccentric 
effects  transferred  to  the  spatial  domain.  This  allows  the  double  integration  of  equa¬ 
tion  5.20c  to  be  accomplished  across  u  and  v  in  the  same  manner  as  the  circular  aperture, 
yielding  f(p'). 


f  (£)  d) ellipse 


f  (P  ) ellipse 


f  {P  ) ellipse 


1 

2n 

1 

2n 


Umax  /  Umax 


r2n 


e2in(u^+vr,)dudv 


z2i-KP'(r  cos  9)  r  dQdr 


'0  Jo 
\ 

J0(2Tipr)  r  dr 


(5.21a) 

(5.21b) 

(5.21c) 


The  normalized  intensity  map  may  be  constructed  by  noting  the  maximum  value  of  Ja  is 
still  1.  Performing  the  normalization  and  squaring  yields 


Iq  (p  )  ellipse 


Io(P  )  ellipse 


r  £\  “I  2 

f0  2  J0{2tt p' r)rdr 

Rx 

fo2  rdr 
'2Ji(ZW)l2 
Dxnp1 


(5.22a) 

(5.22b) 


By  direct  substitution  of  equation  5.19b  the  intensity  map  for  an  elliptic  aperture  becomes 


Io(£)  ??)  ellipse  — 


2Ji  (dx^\I{102  +  {vY 


D^\/  (hC)2  +  (t/)2 


(5.23) 


This  is  a  far  more  efficient  solution  to  implement5,  figure  5.9. 


sFigure  5.9  took  less  than  2  seconds  to  generate  with  the  same  workstation  over  the  same  discretized 
plane  as  figure  5.7. 
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Figure  5.9  The  analytic  solution  to  the  elliptical  aperture  may  easily  be  implemented 
over  a  discretized  image  plane. 


Note  that  equation  5.23  yields  the  circular  aperture  intensity  map  if  a  =  b,  or  along  the 
£  =  0  axis.  Setting  rj  =  0  shows  that  the  aperture  is  essentially  contracted  by  the  ratio  of 
^  in  the  u  direction,  causing  an  expansion  of  the  Bessel  function  along  the  £  axis. 


5.3  Summary 

Analytic  expressions  are  derived  here  for  the  intensity  maps  generated  by  filled  aper¬ 
tures  of  regular  shapes.  Integrating  the  PSF  over  an  area  yields  the  well-known  solutions 
for  square  and  circular  apertures,  but  neither  of  these  known  shapes  conform  to  the  el¬ 
liptical  loci  of  receivers  in  a  satellite  formation,  and  is  therefore  not  quite  suitable  for  use 
as  the  reference  aperture  in  an  interferometry  metric.  A  filled  elliptical  aperture  has  an 
analytic  I0(£,f?),  and  is  derived  over  a  square  image  plane  by  a  conversion  to  eccentric 
polar  coordinates  in  the  aperture  (an  approach  generalized  in  chapter  VII). 
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VI.  Optimization  of  Formation  Geometries 

With  the  metrics  of  chapter  IV  in  hand,  it  is  possible  to  begin  the  task  of  actually  search¬ 
ing  for  optimal  geometries  for  conceptual  systems.  It  is  immediately  apparent  that  such 
searches  are  complex,  even  for  simple  systems.  This  chapter  begins  by  exploring  the 
non-convex  nature  of  the  comparison  metrics,  even  for  the  simplest  formations.  A  one- 
dinrensional  formation  is  solved  explicitly,  and  demonstrates  the  need  for  a  global  search 
routine.  A  genetic  algorithm  is  introduced  as  a  powerful  means  to  find  an  acceptably 
close  approximation  of  the  optimum  receiver  geometry  for  more  complex  formations.  The 
algorithm  is  verified  with  a  simple  idealized  array  set  up  on  the  North  Pole,  then  for  both 
a  constrained  collinear  planar  array  and  a  general  planar  array,  both  with  ideal  viewing 
geometries. 


6.1  Non- Convexity 

The  receiver  placement  problem  as  posed  is  non-convex.  This  may  be  seen  in  the 
simplest  of  all  systems,  a  one-dinrensional  linear  array1. 


Figure  6.1  A  line  of  co-orbiting  satellites  is  the  simplest  possible  formation,  with  p  =  0 
for  each  element  and  b  0  for  all  but  the  first  receiver. 


For  a  one-dinrensional  array  the  intensity  map,  I(£),  simplifies  equation  4.9  to 

j  gl  [I(Ci)  -!<.(&)]  2 

m 


1This  one-dimensional  problem  actually  describes  a  satellite  formation  system,  one  in  which  all  receivers 
have  identical  orbital  parameters  except  for  a  slight  variance  in  mean  anomaly  (an  unsophisticated,  yet 
valid  satellite  formation). 
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From  figure  4.2,  the  distance  chosen  between  each  pair  of  receivers  determines  the  fre¬ 
quency  of  the  (cos)2  function  they  generate.  Adjusting  this  distance  varies  the  cumulative 
difference  between  the  filled  aperture  response  and  the  discrete  receiver  response.  This 
cumulative  error  is  easily  plotted  for  a  single  pair  of  receivers  over  a  gap  varying  from 
[0.0  ...  1.0]  x  D\,  figure  6.2.  Minimizing  this  error  obviously  requires  finding  the  global 


Single  Baseline  Error  Topology 


Figure  6.2  Varying  the  gap  between  discrete  receivers  up  to  a  maximum  of  D\  shows 
two  local  minima,  one  is  obviously  the  global. 


minimum  at  0.33  x  D\  and  avoiding  the  local  minimum  at  0.83  x  D\. 

Placing  a  third  receiver  in  the  aperture  immediately  complicates  the  problem.  By 
locking  one  at  the  origin2,  the  problem  has  one  degree  of  freedom  for  each  of  the  two 


^Fixing  the  position  of  one  receiver  is  actually  an  important  point.  Since  the  discrete  aperture  response 
is  a  function  only  of  the  difference  between  receiver  locations  [35],  allowing  the  full  degrees  of  freedom  for 
all  receivers  creates  a  solution  space  that  resembles  a  Dirichlet  [43]  function.  The  true  minimum  relative 
geometry  exists  in  an  infinite  number  of  positions  within  the  aperture.  This  distinction,  though  discernable 
in  his  results,  was  not  specifically  cited  in  Kong’s  work  [21]. 
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remaining  receivers. 


Error2 


(i-i„)2 


^e— 2i7r(Aufc£) 


sin(7rE>Ag) 

ttDx£, 


2 


Sliding  the  two  receivers  around  from  [0.0  ...  1.0]  x  D\,  summing  the  contributions  of  the 
three  baselines,  and  differencing  them  from  a  filled  one-dimensional  aperture  results  in  a 
two-dimensional  surface  of  error  flux,  with  four  minima,  figure  6.3.  By  symmetry,  two  of 
them  are  copies  of  the  global  at  D\  x  (0.28,0.73)  and  D\  x  (0.73,0.28).  (This  illustration 


Single  Baseline  Permuational  Error  Topology  for  Two  Receiver  Placement 


1 


Figure  6.3  Varying  the  placement  of  two  discrete  receivers  in  a  one-dinrensional  aper¬ 
ture  from  the  origin  out  to  a  maximum  of  D\  adjusts  three  baselines.  The 
aggregate  aperture’s  error  shows  four  local  minima,  two  are  symmetric  copies 
of  the  global. 
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also  highlights  the  limitations  of  graphically  depicting  error  surfaces.  Two  receivers  moving 
in  one-dinrension  use  up  all  available  plotting  dimensions.) 
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6.2  A  Genetic  Algorithm 

Navigating  around  such  multi-dimensional  surfaces  to  approximate  a  global  mini¬ 
mum  immediately  suggests  a  numeric  approach.  Due  to  its  non-convex  nature,  simple 
hill-climbing  (descending)  searches  may  be  immediately  discounted  since  these  approaches 
are  susceptible  to  getting  caught  in  local  minima.  Kong,  Miller,  and  Sedwick  employed  a 
Simulated  Annealing  technique  [21]  that  used  a  global  search,  slowly  reduced  in  scope  by 
lowering  a  theoretical  “temperature”  of  the  system.  Another  popular  technique  (and  the 
one  selected  here)  is  the  use  of  genetic  algorithms.  With  fine  enough  gene  definition,  the 
standard  genetic  process  builds  chromosome  structures  that  represent  the  defining  param¬ 
eters  for  each  receiver’s  paths.  Finite-sized  populations  of  such  structures  accumulate  the 
configurations  yielding  minimal  errors.  Successive  generations  of  crossovers  and  mutations 
vary  these  parameters  within  the  solution  space  and  reevaluate  the  fitness  of  offspring 
members. 

A  conservatively  constructed  set  of  genes  yields  complete  numerical  freedom  for  each 
parameter  to  desired  precision.  This  research  employed  a  simple  real  valued  numeric  and 
scaling  set. 


Gene 

Value 

Gene 

Value 

genel 

1 

gene8 

8 

gene2 

2 

gene9 

9 

gene3 

3 

gene 10 

0.1 

gene4 

4 

gene 11 

0.01 

gene5 

5 

gene 12 

0.001 

gene6 

6 

gene 13 

0.0001 

gene7 

7 

gene 14 

0.00001 

Table  6.1  This  simple  set  of  genes  certainly  spans  the  solution  space,  and  provides  for 
multiple  means  of  defining  each  normalized  chromosome. 

By  constructing  chromosomes  as  sums  and  differences  of  genes  1  through  9,  each  scaled  by 
genes  10  through  14,  a  population  of  coefficients  densely  populating  the  range  [—1.0  ...  1.0] 
may  be  generated  to  an  accuracy  of  five  decimal  places3.  Such  a  population  of  generic 


3 The  range  [—1.0  ...  1.0]  offers  more  flexibility  than  [0.0  ...  1.0]  in  numerically  selecting  parameters. 
The  selection  routine  set  up  to  use  this  range  suffers  less  infant  mortality  from  chromosomes  violating 
constraints. 
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chromosomes  may  be  used  to  scale  each  degree  of  freedom  in  a  formation,  regardless  of  its 
dimensionality.  This  set  is  actually  quite  flexible,  since  the  use  of  addition  and  subtraction 
allows  a  numeric  value  to  be  built  several  redundant  ways,  i.e.  0.98  =  gene9  x  genelO  + 
gene8  x  genell  or  0.98  =  genel  —  gene2  x  genell.  This  flexibility  allows  a  desirable  set 
of  chromosomes  to  be  constructed  in  many  ways  by  the  searching  algorithm.  Just  as  the 
simulated  annealing  technique  of  Kong,  Miller,  and  Sedwick  used  a  heuristically  chosen 
cooling  rate  for  the  system  “temperature”  to  eventually  halt  the  search,  genetic  searches 
may  be  performed  for  an  empirically  chosen  large  number  of  generations  to  arrive  at  an 
acceptably  close  approximation  of  the  optimal  set  of  parameters.  The  definition  of  “large 


Convergence 


Figure  6.4  An  appropriate  number  of  generations  for  performing  the  genetic  algorithm 
may  be  determined  empirically  by  monitoring  the  best  objective  value  of  the 
populations.  By  construction,  each  search  exhibits  a  monotonically  decreas¬ 
ing  objective,  asymptotically  approaching  that  of  the  true  solution. 

number  of  generations”  grows  with  the  degrees  of  freedom  in  the  formation,  but  may  be 
determined  by  monitoring  the  best  fitness  value  for  the  populations  during  the  course  of 
the  search,  figure  6.4.  As  the  search  progresses,  by  construction  the  best  fitness  value 
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exhibits  a  monotonically  decreasing  behavior.  With  each  generation  of  the  search,  the 
algorithm  refines  the  geometry  of  the  formation  based  on  this  optimization  process.  Since 
this  algorithm  is  by  nature  a  stochastic  search,  a  truly  exhaustive  search  would  require 
performing  the  search  multiple  times  and  performing  statistical  validation  of  the  results. 
However,  this  research  only  seeks  validation  of  the  problem  formulation  and  the  practicality 
of  approximated  solutions.  Note  that  the  solution  surface  of  Jimage  remains  analytic  for 
any  size  formation,  figure  6.5.  Though  “ugly,”  the  summation  formulation  for  Ip (£,77,  t) 
has  no  discontinuities.  Likewise,  derivatives  for  the  Bessel  function  formulation  of  I0(£,  77) 
all  exist4.  Therefore  a  single  search  of  each  smooth  error  surface  is  performed  to  find  an 


Elliptical  Mean  Permuational  Error  Topology  for  Two  Receiver  Placement 


Figure  6.5  A  similar  technique  to  that  used  in  figure  6.3  may  be  used  to  generate  a 
surface  plot  of  Jimage  for  a  two  degree  of  freedom  formation.  The  loci  and 
performance  of  a  planar  collinear  formation  (section  6.4)  with  three  receivers 
is  determined  by  the  selection  of  pH  for  the  two  movable  vehicles. 

acceptable  approximate  solution  among  the  many  non-unique  possible  solutions. 

4 Any  standard  mathematics  handbook  [10]  gives  the  expanded  approximation  of  the  Bessel  function  in 
^  ( _ i)fcf  ^'jrl+2fc 

a  form  like  J(x)n  =  fc!r(w+fc+i)  ’  which  is  another  exercise  in  differentiation,  but  the  derivatives 

exist. 
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The  overall  process  for  this  genetic  algorithm  is  summarized  in  the  logic  flow  of 


table  6.2.  (A  more  extensive  description  is  included  as  Appendix  E.) 


Step  Operation 

1  Define  gene  pool. 

2  Define  initial  chromosomes. 

3  Generate  population  of  coefficient  chromosomes  for  each  degree  of  freedom. 

4  Evaluate  fitness  (normalized  error)  of  initial  populations,  equation  4.18. 

5  If  population  is  too  large  trim  chromosomes  with  highest  fitness. 

6  Select  chromosomes  for  genetic  operation,  random,  weighted 
for  chromosomes  with  lower  fitness  values. 

7  Perform  crossover/mutation  operations  on  chromosomes. 

8  Convert  new  chromosomes  to  numerical  values. 

9  If  constraints  satisfied,  calculate  fitness  of  new  chromosomes, 
otherwise  delete. 

10  If  less  than  maximum  number  of  generations,  go  to  5. 

11  Manual  check  for  convergence. 

(yes,  continue;  no,  increase  maximum  generations,  go  to  5) 

12  Generate  array  performance  outputs. 

Table  6.2  This  algorithm  summarizes  the  standard  programming  logic  [34]  used  to  per¬ 
form  the  genetic  searches  in  this  investigation. 


The  validity  of  this  single  search  process  is  supported  by  an  example  stochastic  search. 
Figure  6.4  is  actually  one  of  ten  repeat  searches  performed  for  a  five-receiver  general  planar 
formation  (introduced  in  section  6.5),  each  with  different  sets  of  initial  conditions.  The 
results  of  the  searches  showed  very  little  dispersion  about  the  mean  (standard  deviation  of 
cr  =  0.0000520).  This  small  dispersion  is  illustrated  by  the  tiny  2er  error  bar  drawn  above 


Search 

J  image 

Search 

J  image 

1 

0.00046534 

6 

0.00049152 

2 

0.00057438 

7 

0.00054121 

3 

0.00052650 

8 

0.00059671 

4 

0.00058134 

9 

0.00051591 

5 

0.00045729 

10 

0.00059498 

a 

0.00005200 

Table  6.3  Actually  performing  multiple  searches  for  the  smallest  formation  in  section  6.5 
shows  a  very  narrow  distribution  about  the  mean. 

the  final  generation  in  figure  6.4.  This  interval  statistically  encompasses  95.449%  of  the 
numerically  approximated  least  upper  bounds  [48]  on  the  true  minimum  Jimage- 
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Since  the  error  surface  of  Jimage  is  analytic  for  any  size  formation  and  the  example  stochastic 
search  for  a  small  formation  demonstrates  a  very  narrow  dispersion  of  solutions,  performing 
single  searches  for  larger  formations  should  not  risk  drastic  deviations  away  from  the  true 
minimum  aggregate  error,  and  is  deemed  adequate  for  valid  solution  approximations. 

This  research  used  a  modified  form  of  the  program  architecture  set  forth  by  the 
Genetic  Search  Toolbox®  by  Optimal  Synthesis  Inc.®  [26].  The  toolbox  provides  functions 
for  generating  genes,  chromosomes,  and  populations  as  character  string  representations. 
Genetic  operations  (crossovers,  mutations)  are  performed  on  the  character  strings,  which 
must  then  be  evaluated  to  yield  numerical  values  for  the  fitness  (aggregate  error).  The  basic 
search  structure  recommended  by  the  toolbox  authors  was  expanded  to  provide  progress 
updates,  a  graphical  convergence  display,  and  user-controllable  extension  of  the  search 
generations  seeking  the  behavior  described  by  Rahmat-Samii  and  Michielssen  [34].  Each 
search  is  performed  until  successive  reductions  in  the  objective  function  appear  to  be  less 
than  10%  of  the  latest  Jimage- 

6.3  A  Simplified  Polar  Array  Example 

The  basic  problem  definition  and  solution  technique  may  be  illustrated  through  anal¬ 
ysis  of  an  interferometric  array  with  a  simple  geometry.  A  terrestrial  array  distributed 


Figure  6.6  An  array  of  receivers  distributed  about  a  pole  offers  very  simple  propagation 
and  path  loci. 

about  the  North  Pole  offers  perfectly  circular  receiver  paths  and  elementary  propagation, 
figure  6.6.  For  imaging  high  declination  targets,  a  dwell  time  of  half  a  sidereal  day  may 
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be  assumed,  which  further  reduces  the  dimensionality  of  the  problem.  Since  a  rotation 
of  180°  swings  each  baseline  through  all  possible  orientations,  a  linear  array,  determined 
only  by  radial  distance,  eventually  offers  all  the  performance  of  more  complex  geometries 
(though  instantaneous  or  “snapshot”  performance  will  be  poor).  The  circular  paths  also 
suggest  the  simplification  of  comparing  the  array  to  a  filled  circular  aperture,  with  the 
known  solution  from  section  5.2.2. 

Consider  setting  up  such  an  array  of  five  receivers,  one  at  the  pole  and  the  others  ar¬ 
rayed  about  it.  This  system  has  four  degrees  of  freedom,  requiring  four  generic  populations, 
each  selecting  a  radial  distance,  [—1.0 ...  1.0]  x  For  purposes  of  non-dimensionality, 
select  the  maximum  array  diameter5  as  D\  =  100  wavelengths.  Monitoring  the  best  ap- 


Convergence  10  Baseline  Lengths  Propagated  Best  Positions 
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Figure  6.7  After  20,000  generations,  a  solution  for  five  receivers  tracing  out  paths  about 
the  North  Pole  minimizes  its  aggregate  error  as  compared  to  a  circular  aper¬ 
tures.  The  baseline  distribution  demonstrates  the  natural  inclusion  of  near¬ 
maximum  length  baselines,  with  enough  short  baselines  to  reduce  sidelobe 
development  of  the  high  spatial  frequencies. 

proxinration  in  the  populations  at  the  end  of  each  generation  shows  that  a  single  search  of 
20,000  generations  appears  sufficient  to  achieve  a  near-optimal  solution  for  receiver  spac- 
ings  whose  performance  closely  mimic  the  circular  aperture.  The  lower  right  graphic  in 

5  For  purposes  of  consistency,  the  convention  of  setting  the  largest  dimension  of  an  aperture  to 
D\  =  100  wavelengths  is  followed  throughout  this  work  whenever  an  arbitrary  aperture  size  must  be 
selected. 
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figure  6.7  shows  the  performance  of  such  an  array  over  a  180°  arc.  Inspection  of  the  cho¬ 
sen  radii  for  the  array  show  an  almost  perfect  distribution  of  baseline  lengths,  including 
near-maximum  length  (for  high  spatial  frequencies)  and  one  of  minimal  length  (necessary 
to  minimize  the  sidelobes  in  the  composite  intensity  map).  Even  for  this  simplified  prob¬ 
lem,  this  solution  exists  in  several  locations  on  the  J image  surface.  Simply  trading  the 
assignment  of  p  for  pairs  of  receivers  gives  several  equivalent  solutions. 

Four  receivers  is  not  quite  sufficient  to  demonstrate  receiver  placement  trends,  so  the 
array  was  expanded  to  nine  receivers,  figure  6.8.  This  creates  an  eight  degree  of  freedom 
problem,  and  the  solution  technique  is  the  same,  it  just  requires  50,000  generations  for  con- 
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Figure  6.8  After  50,000  generations,  an  optimized  array  of  nine  receivers  about  the  North 
Pole  takes  advantage  of  its  ability  to  more  densely  fill  the  {u  v)  plane. 

vergence.  (Recall  that  at  this  complexity,  doubling  the  degrees  of  freedom  creates  nearly 
four  times  the  number  of  baselines  of  the  smaller  array.  This  not  only  increases  the  run 
time  for  the  simulation  to  propagate  the  receivers  over  180°,  but  vastly  increases  the  num¬ 
ber  of  possible  relative  receiver  locations.  Only  judicious  selection  of  the  correct  members 
of  each  population  for  crossover/mutation  or  removal  keeps  the  search  to  a  practical  dura- 
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tion6.  The  optimized  distribution  shows  the  necessary  inclusion  of  many  maximum-length 
baselines,  but  over  half  of  the  permutation  baselines  are  still  chosen  to  yield  the  shorter 
lengths,  giving  the  low  spatial  frequencies  needed  to  generate  a  complete  image  (review 
figure  3.9). 

6.4  Optimizing  a  Planar  Collinear  Formation 

The  notion  from  section  4.4  of  optimizing  receiver  motion  over  time  necessitates 
including  path  shape  information  into  the  selection  of  a  reference  aperture.  The  simple 
circular  paths  of  the  previous  section  rarely  happen  for  actual  systems.  From  the  descrip¬ 
tions  of  Wiesel  [54],  and  inspection  of  equations  2.2,  receivers  in  a  formation  of  satellites 
trace  out  ellipses  in  the  (ex,ey)  plane  of  their  local  coordinate  system.  For  high  relative 
declination  targets,  this  defines  the  (u  v)  plane  for  all  practical  purposes  (an  assumption 
that  is  relaxed  in  the  next  chapter).  Therefore  for  any  choice  of  rnH  and  nH  defining 
the  plane  for  the  formation  and  for  any  choice  of  the  two  remaining  Hill’s  parameters  pH 
and  0H,  the  receivers’  paths  trace  out  ellipses  with  a  2:1  ratio  of  major  to  minor  axes 
(eccentricity  of  ^). 


Dx 


Figure  6.9  A  planar  formation  may  be  further  simplified  by  assuming  a  common  0H for 
each  receiver.  A  180°  loci  rotates  these  baselines  through  every  possible 
orientation. 

6The  optimization  depicted  in  figure  6.8  took  nearly  five  days  to  converge  on  a  1.7  GHz  workstation. 
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Choosing  a  reference  aperture  that  mimics  these  loci  avoids  the  inconsistencies  de¬ 
scribed  in  section  4.3.1.  An  initial  problem  may  be  simplified  by  choosing  a  comparison 
time  of  half  an  orbit  [32],  figure  6.9.  This  removes  a  degree  of  freedom  from  each  receiver 
(just  as  the  half-day  period  did  for  the  polar  array)  since  each  baseline  rotates  through  all 
possible  orientations.  Such  a  linear  array’s  receiver  paths  may  be  specified  by  a  selection 
of  pH  only.  Again,  this  simplification  would  never  be  used  in  a  practical  system,  since  the 
snapshot  performance  is  unreasonably  poor,  figure  6.10.  (Though  the  baselines’  lengths 


Snapshot  Performance,  4  Receivers 


eta  (rad) 

Figure  6.10  Simplifying  a  planar  array  to  stay  in  a  line  within  the  (ex,ey)  plane  may 
yield  adequate  performance  over  some  dwell  time,  but  the  instantaneous  or 
“snapshot”  performance  is  only  a  one-dinrensional  distribution. 

also  change  over  the  course  of  the  dwell  time,  each  length/orientation  will  eventually  be 
sampled  in  a  180°  arc.) 

Beginning  with  a  modest  formation  of  nine  receivers,  the  problem  has  eight  degrees 
of  freedom,  eight  values  of  pH  chosen  from  the  range7  [—1.0...  1.0]  x  Figure  6.11 

7Recall  the  definition  of  pH  from  Hills  equations. 
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shows  just  such  an  optimized  formation.  Note  the  distribution  of  baseline  lengths  in  this 


Convergence  Distribution  of  36  Baseline  Lengths 


Figure  6.11  Optimizing  a  nine  receiver  formation  for  half  an  orbit  puts  the  dimensional¬ 
ity  of  the  problem  to  just  eight,  and  40,000  generations  provide  convergence. 


figure.  This  is  perhaps  the  most  informative  formation  description,  since  it  determines 
the  sampled  spatial  frequencies  at  each  discrete  time.  As  expected,  the  solution  locations 
include  provisions  for  several  maximum  length  baselines,  and  a  nearly  uniform  sampling 
of  the  other  lengths  over  the  (u  v)  plane. 

To  see  the  utility  of  comparing  the  formation’s  performance  to  a  filled  elliptical 
reference,  temporarily  ignore  the  shape  of  the  receiver  paths,  and  attempt  to  drive  the 
formation  to  mimic  a  simple  circular  aperture.  The  choice  for  such  an  aperture  falls 
to  either  a  smaller,  inscribed  circular  aperture,  or  a  larger,  circumscribed  aperture  (see 
figure  4.3).  Figure  6.12  shows  the  same  nine-receiver  formation  just  created,  but  optimized 
as  compared  to  a  small  circular  aperture.  (Note  the  smaller  reference  aperture  generates  a 
wider  peak  in  the  image  plane.)  Hampered  by  a  poor  reference  standard,  the  optimization 
accepts  suboptimal  geometries  that  generate  the  wider  peak  of  its  reference,  apparent  in 
the  reduced  number  of  near-maximum  length  permutational  baselines  generated. 

Comparing  the  distributed  formation  to  a  larger,  circumscribed  circular  aperture 
fares  better  in  the  optimization  as  compared  to  the  inscribed  circular  reference.  Parts 
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Convergence  Distribution  of  36  Baseline  Lengths 
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Figure  6.12  Optimizing  the  nine  receiver  formation  compared  to  a  smaller  inscribed 
circular  aperture  sized  to  match  the  minor  axis  yields  a  formation  that  fails 
to  optimally  construct  baselines  filling  the  (u  v)  plane. 


of  the  optimal  paths  in  the  reference  are  not  reachable  by  the  formation  members  but  a 
similar  solution  as  the  smaller  reference  is  eventually  found,  figure  6.13. 
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Figure  6.13  Optimizing  the  same  nine  receiver  formation  compared  to  a  larger  circular 
aperture  sized  to  match  the  major  axis  eventually  yields  a  comparable  result 
to  the  elliptical  reference. 
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This  comparison  suggests  reduction  to  a  lower  ordered  problem  to  gain  some  insight 
into  what  the  optimization  is  generating.  Keeping  the  relative  linear  geometry  at  each 
time  step,  but  halving  the  number  of  non-fixed  receivers  decreases  the  degrees  of  freedom 
to  four  (for  a  five- receiver  formation),  figure  6.14.  With  the  reduced  number  of  baselines, 
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Figure  6.14  Optimizing  a  five  receiver  formation  offers  a  lower-dimensionality  compar¬ 
ison  of  comparison  metrics.  This  optimization  required  20,000  generations 
for  convergence,  and  displays  the  same,  though  sparser  distribution  as  fig¬ 
ure  6.11. 


the  overall  distribution  of  lengths  include  provisions  for  a  maximum-length  sample  and  a 
nearly  uniform  sampling  of  the  rest  over  the  (u  v)  plane. 

Driving  the  smaller  formation  to  mimic  a  small  circular  aperture’s  performance  yields 
a  similar  distribution.  Figure  6.15  shows  the  same  five-receiver  formation  just  created 
but  optimized  as  compared  to  an  inscribed  circular  aperture.  The  wider  central  peak 
generated  by  a  smaller  aperture  again  allows  the  optimization  search  to  more  readily 
accept  suboptimal  geometries,  apparent  in  the  lack  of  a  maximum-length  permutational 
baseline. 

Performing  the  same  optimization  with  a  large  circular  aperture  has  similar  results  to 
the  elliptical  reference,  figure  6.16.  In  fact,  at  the  lower  dimensionality  it  is  possible  to  see 
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Figure  6.15 


Figure  6.16 
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Optimizing  the  five  receiver  formation  referenced  to  a  smaller  circular  aper¬ 
ture  still  allows  it  to  accept  suboptimal  performance  due  to  the  wider  central 
peak  of  a  small  aperture. 
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Optimizing  the  same  five  receiver  formation  compared  to  a  larger  circular 
aperture  sized  to  match  the  major  axis  eventually  yields  a  comparable  result 
to  the  elliptical  reference. 
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the  nearly  perfect  one-to-one  correspondence  with  the  solution  offered  by  the  elliptically 
referenced  solution.  The  optimization  based  on  the  larger  circular  reference  eventually 
drives  the  receiver  placement  out  to  the  constraints  of  the  dynamics.  Note,  however,  that 
both  figures  6.15  and  6.16  show  a  much  more  sluggish  convergence  performance  than  the 
optimization  based  on  the  more  appropriate  elliptical  reference. 

6.5  Optimizing  a  General  Planar  Formation 

Though  the  linear  formations  of  section  6.4  were  instructive  for  simple  solutions  and 
comparisons,  they  suffer  from  poor  “snapshot”  performance.  Since  each  would  require  very 
long  dwell  times  to  yield  acceptable  images,  it  is  unlikely  such  a  geometry  would  ever  be 
selected.  A  more  realistic  formation  allows  complete  freedom  of  receiver  placement  within 
the  formation  plane[30],  figure  6.17. 

IT 
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Figure  6.17  A  planar  formation  provides  a  natural  definition  of  the  ( u  v )  plane.  Re¬ 
ceivers  in  the  formation  trace  out  simple  2:1  ellipses,  but  the  permutational 
baselines  may  have  fairly  varied  orientations,  depending  on  the  relative  lo¬ 
cations  of  the  receivers. 
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Consider  another  five  receiver  formation,  figure  6.18,  with  one  receiver  locked  at  the 
origin.  Each  receiver  may  be  placed  by  selecting  both  its  pH  and  0H  at  the  expense  of 
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Figure  6.18  Optimizing  the  five  receiver  formation  of  section  6.4  may  be  readily  en¬ 
hanced  to  allow  two  degrees  of  freedom  per  receiver.  Placing  one  receiver  at 
the  origin  leaves  eight  degrees  of  freedom  in  the  formation.  The  enhanced 
performance  over  summing  the  linear  arrays  is  immediately  apparent. 

doubling  the  degrees  of  freedom.  The  formation’s  instantaneous  performance  is  now  quite 
reasonable.  The  simulation  time  is  also  reasonable.  Though  eight  populations  must  be 
generated,  there  are  still  only  four  receivers  to  propagate  and  ten  baselines  to  be  generated 
for  the  aggregate  performance  at  each  time  step.  Since  an  initial  0H  may  be  chosen  for 
each  receiver,  the  dwell  time  may  be  halved  to  investigate  which  baseline  orientations  are 
the  most  beneficial. 

If  the  genetic  algorithm  is  allowed  to  select  pH  =  [—1.0...  1.0]  x  ^  and  0H  = 
[—1.0. . .  1.0]  x  90°,  a  dwell  time  of  one  quarter  orbit  offers  all  possible  orientations,  but 
not  all  possible  baseline  lengths  (a  more  realistic  limitation).  Figure  6.18  depicts  such  an 
optimized  formation.  Since  simple  initial  baseline  length  is  no  longer  sufficient  to  illustrate 
the  sampled  spatial  frequencies,  figure  6.18’s  distribution  diagram  has  been  modified  to  a 
polar  depiction  to  show  the  length  and  orientation  of  all  initial  permutational  baselines. 
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The  inclusion  of  near-maximum  length  baselines  is  apparent,  but  the  shorter  comparison 
time  drives  them  to  be  repeated  at  several  orientations. 


Figures  6.19  through  6.22  show  the  successive  improvements  in  formation  perfor¬ 
mance  by  simply  adding  receivers.  As  the  formation  grows,  the  aggregate  aperture  ap- 
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Figure  6.19  A  six-receiver  formation  has  ten  degrees  of  freedom.  After  40,000  genera¬ 
tions,  the  optimized  formation  has  essentially  perfect  performance. 


proximates  the  filled  reference  aperture  very  closely  and  this  behavior  is  quantified  in 
the  tiny  values  of  the  objective  function.  Expanding  the  formation  out  to  7,  9,  and  11 
receivers  yields  only  marginal  improvements  in  the  near-perfect  performance,  but  does 
require  extended  numbers  of  generations  for  convergence.  This  outstanding  performance 
will  of  course  be  degraded  by  reducing  the  target  dwell  time  to  less  than  the  idealized  90° 
of  orbit  arc. 
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Figure  6.20 


A  seven-receiver  formation 


has  twelve  degrees  of  freedom. 
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Figure  6.21 


A  nine-receiver  formation  has  sixteen 


degrees  of  freedom. 
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Figure  6.22  An  eleven-receiver  formation  has  twenty  degrees  of  freedom  (and  requires 
100,000  generations  for  adequate  convergence).  Its  overall  performance  over 
a  90°  sweep  is  comparable  to  smaller  formations. 


Actual  computing  time  for  these  larger  numbers  of  generations  is  actually  only 
slightly  longer  over  smaller  formations  due  to  the  increased  infant  mortality  in  new  gen¬ 
erations.  With  larger  numbers  of  populations,  the  early  generations  have  an  increased 
probability  that  one  of  its  chromosomes  will  violate  the  size  constraint  and  be  eliminated. 
This  skips  initiating  a  simulation  for  that  prospective  member. 

Overall  the  benefits  of  increasing  the  formation  size  results  in  a  diminishing  error 
compared  to  the  filled  aperture,  figure  6.23.  Appendix  F  has  the  complete  results  of 
optimized  formations  of  three  and  four  receivers. 


6.6  Summary 

Refining  a  formation  based  on  an  optimization  process  using  the  metrics  of  chapter  IV 
is  a  non-convex  problem.  The  solutions  available  are  non-unique,  and  may  be  determined 
by  employing  the  global  search  techniques  of  genetic  algorithms.  The  numerically  derived 
solutions  offer  performances  that  show  improvement  with  both  increased  reference  aperture 
fidelity  and  increased  formation  density.  These  initial  problems  were  posed  with  reduced 
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Figure  6.23  Increasing  the  number  of  receivers  in  the  formation  drives  the  performance 
closer  to  that  of  the  filled  aperture,  with  predictably  diminishing  improve¬ 
ments  with  each  new  receiver. 


complexity  and  dimensionality  by  assuming  ideal  viewing  geometries  (the  targets  were 
assumed  to  be  directly  normal  to  the  formation  planes) .  This  assumption  is  finally  relaxed 
in  the  next  chapter. 
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VII.  Optimizing  for  Inertial  Targets 

With  the  idealized  viewing  geometry  constructions  of  the  previous  chapters  in  hand,  it  is 
possible  to  extend  the  solution  search  to  more  realistic  premises:  practical  formations  and 
real  targets[31].  This  chapter  defines  a  notional  reference  orbit  for  a  formation  and  locates 
four  example  stellar  targets.  Both  the  aggregate  aperture  loci  and  reference  aperture 
shapes  are  then  warped  to  model  the  distortions  caused  by  the  non-ideal  viewing  geometry 
from  each  target.  A  growing  series  of  formations  is  then  re-solved  for  each  target  using  the 
genetic  algorithm  of  chapter  VI. 

7. 1  Inertial  Targets 

Consider  a  low-flying,  high  inclination  formation,  one  that  offers  cheap  launch  costs 
and  excellent  viewing  geometries  for  objects  in  the  vicinity  of  the  celestial  equator  (and 
the  orbit  normal  vector)1. 


Figure  7.1  A  low  altitude,  high  inclination  formation  offers  excellent  viewing  geometries 
about  the  equatorial  regions. 


1  Using  using  short  enough  wavelengths,  such  a  formation  might  also  be  used  for  viewing  high  altitude, 
low  inclination  Earth-orbiting  objects. 
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For  consistency,  select  a  simple  low  altitude  circular  orbit  inclined  at  1  radian,  and  an 
ascending  node  at  17  =  150°.  Such  an  orbit  offers  relatively  simple  propagation  and  good 
viewing  geometry  for  objects  within  90°  of  the  Celestial  Origin.  Choose  to  start  the 
imaging  arc  at  an  arbitrary  true  (mean)  anomaly  of  v  =  40°.  The  main  complication  in 
optimizing  a  formation  about  such  a  reference  orbit  for  an  actual  inertial  target  is  that 
the  assumption  of  near  90°  relative  declination  is  no  longer  valid.  The  {u  v  w)  coordinate 
frame  is  determined  by  the  target  with  the  w  axis  exactly  aligned  with  the  target  look 
vector.  The  local  {x  y  z)  frame  is  determined  by  the  formation  dynamics.  Note  that  the 
two  frames  are  not  aligned.  To  illustrate  the  technique  required  to  align  the  reference 
frames,  it  is  helpful  to  have  a  few  example  targets. 

Celestial  objects  are  located  by  coordinates  of  declination  (5,  and  either  hour  angle, 
H ,  or  right  ascension,  $.  Since  $  is  measured  eastward  from  the  same  reference  as  the 
right  ascension  of  the  ascending  node,  Q,  it  is  the  logical  choice  here  over  H.  Four  celestial 
targets  [47],  chosen  for  proximity  to  the  Celestial  Equator  and  brightness  are  given  in 
table  7.1.  The  chosen  orbit  offers  good  viewing  geometry  for  the  “nearby”  Andromeda 
galaxy  and  the  stars  Betelgeuse  and  Rigel,  both  in  the  constellation  of  Orion.  Regulus, 
in  the  constellation  Leo,  is  deliberately  included  in  this  evaluation  for  its  poor  viewing 
geometry. 


Object 

Right  Ascension 

Declination 

M31  (Andromeda) 

10.675° 

41.267° 

HR2061  (Betelgeuse) 

88.793° 

7.407° 

HR1713  (Rigel) 

78.634° 

-8.202° 

HR3982  (Regulus) 

152.093° 

11.967° 

Table  7.1  Each  of  these  bright  celestial  objects  is  located  by  the  coordinates  of  right 
ascension,  $,  and  declination,  5. 


7.2  Orienting  the  Aggregate  Aperture 

To  map  out  the  effective  intensity  map  generated  by  the  aggregate  aperture  of  a 
formation,  the  w  axis  of  the  spatial  frequency  coordinate  system  must  be  aligned  with 
the  direction  of  the  intended  target.  This  carries  the  minor  assumption  from  chapter  IV 
that  the  system  uses  sufficiently  long  wavelengths  that  phase  may  be  measured  and  all 
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receivers’  signals  mapped  to  the  same  (u  v)  plane.  This  alignment  requires  two  rotations 
through  the  angles  6  and  <f>,  but  these  angles  are  in  the  inertial  frame,  not  the  non-inertial 
Clohessy-Wiltshire  frame.  Therefore  the  receiver  locations  in  the  non-inertial  coordinate 
system  must  first  be  mapped  out  to  the  inertial  system. 
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(x  y  z) 
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(u  V  w ) 
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(7.1) 


This  transformation  is  the  reverse  of  that  used  by  Wiesel  [54],  summarized  in  Ap¬ 
pendix  C,  to  map  an  inertial  Hamiltonian  out  to  the  non-inertial  frame.  For  a  distant 
source,  it  is  not  necessary  to  perform  the  origin  translation  by  N0  =  (X'Q ,  Y'0.  Z'Q ,  P'Xo,  Py0-,  P'z0)- 


Figure  7.2  A  stellar  source  located  at  (4>,  5)  determines  the  orientation  of  the  interferom¬ 
eter  coordinate  system  (u  v  w).  The  baselines  between  each  pair  of  receivers 
in  inertial  coordinates  can  be  resolved  into  this  system  by  two  rotations  of 
|  +  <f>  and  |  —  5. 

From  figure  7.2,  this  leaves  the  rotations  of 
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See  Appendix  C  for  general  definitions  of  the  inertial  to  nodal  frame  rotation  matrix, 
TZz  (ignoring  17  for  the  short  duration  dwell  time),  and  the  nodal  to  local  C-W  frame 
rotation  matrix,  TZ.  Setting  up  the  inertial  rotations2  and  inserting  the  rotations  from  the 
Clohessy-Wiltshire  frame  gives 
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Simplifying, 
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This  adjusts  vectors  in  the  local  (x  y  z)  system  to  be  remapped  to  the  ( u  v  w)  axes,  an 
inertial  frame  for  an  inertial  target.  Using  the  simple  geometry  of  figure  7.1,  TZz  and  7 Z 
may  be  constructed  explicitly  in  terms  of  the  Keplerian  elements  of  the  reference  orbit. 
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Simplifying, 


2Thompson  [46]  gives  a  related  formulation  of  these  rotations  to  translate  stellar  coordinates  to  a  ter¬ 
restrial  interferometer  frame.  These  may  be  applied  to  the  relative  coordinates  of  a  formation  provided 
they  are  preceded  by  the  1Z%1ZT  rotations. 
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Aligning  the  local  coordinate  system  with  the  look  vector  to  a  stellar  target  has  the  effect 
of  projecting  the  elliptical  formation  onto  the  plane  orthogonal  to  the  look  vector.  This 
skews  the  perfect  elliptical  loci  of  the  formation  paths,  figure  7.3.  With  the  propagated 
Loci  in  Formation  Frame  Seen  from  Andromeda 


Seen  from  Betelgeuse  Seen  from  Rigel  Seen  from  Regulus 


u  (wavelengths)  u  (wavelengths)  u  (wavelengths) 


Figure  7.3  Moving  the  viewer’s  vantage  point  over  to  the  look  vector  from  each  of  the 
four  stellar  targets  demonstrates  how  the  formation  loci  appear  skewed  when 
viewed  from  a  source  not  directly  above  the  formation. 

receiver  locations  mapped  to  the  inertial  (u  v)  frame,  the  aggregate  aperture  may  be 
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synthesized  using  the  same  permutations  developed  in  equation  4.12  and  intensity  map  of 
equation  4.13. 

7.3  Orienting  the  Reference  Aperture 

The  mapping  of  equation  7.6  gives  the  appropriate  length  and  orientation  for  each 
individual  baseline  used  to  build  the  aggregate  aperture  intensity,  Ip  (£,  7/,f),  but  the  ref¬ 
erence  intensity  map,  1 0(£,  77),  must  also  be  altered  to  reflect  the  true  relative  geometry. 
Otherwise  the  optimization  is  hampered  by  seeking  loci  where  the  formation  dynamics  do 
not  allow  them.  This  alteration  may  be  accomplished  by  generalizing  the  derivation  of 
equation  5.23. 

Consider  the  basic  polar  and  rectangular  representations  of  the  PSF  from  equa¬ 
tions  5.18, 


PSFellipse(r,  9,  p')  =  eW(rcos(*))  (7.7a.) 

PSFelIipse(u,u,C,??)  =  e2i<<+vv)  (7.7b) 

which  allowed  the  intensity  map  symmetric  in  9  to  be  skewed  in  the  spatial  domain  image 
plane  by  defining  the  eccentric  radius,  p generated  by  adjusting  its  components’  rectan¬ 
gular  coordinates 

pi=\f(W^'  <7-8) 

An  elliptical  aperture,  viewed  from  a  direction  off-normal  to  the  relative  orbit  plane 
will  have  some  further  distortion  of  its  apparent  shape.  Generalizing  the  technique  used  in 
section  5.2.3,  this  shape  may  be  determined  in  the  spatial  frequency  domain  and  transferred 
through  the  PSF  to  define  a  warped  eccentric  radius,  p",  in  the  spatial  domain.  Each  of 
the  rectangular  components  of  the  spatial  frequency  radius  may  be  stretched  or  contracted 
individually  by  some  coefficient  to  give  the  warped  eccentric  radius,  r" .  Initially  assume 
they  are  independent  and  define  them  generically  as  (f  in  the  u  direction,  and  if  in  the  v 
direction.  Obviously  each  will  be  some  function  of  the  rotations  from  the  reference  orbit 
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to  the  target  location. 


4>  =  ®,6) 

^  =  f2(i,tt,v,  $,6) 

These  may  be  incorporated  into  the  definition  of  the  warped  eccentric  radius  as  shown 
in  figure  7.4.  This  warped  ellipse  may  be  modelled  as  the  projection  of  each  transformed 


Figure  7.4  The  elliptical  aperture’s  relative  coordinate  system  representation  is  distorted 
when  viewed  from  the  look  angles  (<I>,e>)  to  the  target.  The  distorted  ellipse 
may  be  generated  by  rotations  through  <h,  5,  and  the  reference  orbit  orienta¬ 
tion  angles  u,  i,  and  14. 

radius  onto  the  inertial  (ft  v)  plane.  Using  the  rotations  of  equation  7.6  the  generalized 
PSF eiiipSe*  for  the  image  plane  coordinates  may  be  built. 


F>SFeiiipSe*(ft,  v .  .  tj ) 


PSFeiiipse*(u,  v .  .  // ) 


(7.10a) 


(7.10b) 


For  shorthand,  abbreviate  the  rotation  matrices  from  equation  7.6  as  y  Sx  J ,  y  (|  +  —  14)  ^ 

ix  ,  and  vz  ■  Inserting  these  and  the  contraction  of  ^  into  the  PSF  for  [u  v  w]T 
yields 
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where  the  warped  spatial  coordinates  (£,  iy)*  are  defined  by  transferring  the  rotations  from 
the  spatial  frequency  domain  to  the  spatial  domain. 
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The  PSF  may  now  be  expressed  in  its  polar  form  in  terms  of  the  warped  eccentric  radius, 

PSFellipse*(r,  #,  p")  =  (7.!3) 


p"  as 


This  is  the  same  form  of  the  PSF  as  equation  5.18,  with  the  additional  distortions  all 
contained  in  the  formulation  of  p" ,  warping  the  symmetric  (in  6)  output.  Therefore  p" 
may  be  formally  defined 
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By  setting  the  third  element  of 


l  V  0 


to  zero,  p"  warps  the  known  output  of  a  filled 


L  J  * 

planar  circular  spatial  frequency  aperture,  and  generates  its  projection  onto  the  spatial 
(£,  iy)*  image  plane.  Expanding  the  terms  in  p"  and  using  ||  ■  1 1 2  for  the  vector  magnitude 
allows 
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where  the  subscript  2x2  indicates  just  the  upper  left  four  elements  of  the  rotation  matrix3. 
For  compactness,  define  the  series  of  rotations  as  the  matrix  W 
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Extract  only  the  projection  of  the  image  intensity  map  onto  the  (£,  1))*  plane  by  pulling 
out  the  upper  left  2x2  portion  of  W  to  define  W*. 


_  _ 
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W*(l,  1)  =  Q)  {cos  (f  +  —  12)  cos  z'  +  sin  (|  +  —  12)  cos i sin  z/} 

W*(l,  2)  =  —  cos  (|  +  $  —  12)  sin  z/  +  sin  (|  +  $  —  12)  cos  i cos  v 

>V*(2, 1)  =  Q)  {  —  sin  6  sin  (|  +  $  —  12)  cos  v  +  sin  u  (sin  6  cos  (|  +  —  12)  cosi  +  cos  6 sin  i) } 

W*(2, 2)  =  sin  5 sin  (|  +  —  12)  sin  v  +  cos  v  (sin  <5 cos  (|  +  $  —  12)  cos  i  +  cos  5 sin i) 

(7.17) 


3Since  the  full  3x3  rotation  matrix  is  full  rank,  using  the  full  matrix  multiplied  by  [  £  r/  0  ]  keeps 
the  w  components.  This  does  not  project  the  radius  onto  the  (£,??)*  plane  with  the  resulting  contractions, 
but  preserves  the  original  magnitude.  Using  the  2x2  upper  left  portion  generates  this  projection. 
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Now  equation  5.23  may  finally  be  brought  forward  and  generalized  by  direct  substitution 
of  p"  for  p',  to  give  I0(£,7?)*  as 
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This  describes  the  intensity  map  generated  by  an  aperture  of  elliptical  shape,  following  a 
reference  orbit  with  initial  true  anomaly,  inclination,  and  right  ascension  of  the  ascending 
node  of  (is,  i ,  0) ,  as  viewed  from  an  inertial  source  at  right  ascension  and  declination  (<h,  S). 


7-4  Viewing  Geometry  Effects  on  I0(£,f?) 


Inspection  of  the  upper  left  2  x  2  of  each  rotation  in  W  reveals  the  effects  of  each 
on  (£,77)*.  Rotations  through  is  reorient  the  formation,  but  do  nothing  to  change  its  size, 
since 
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The  formation’s  apparent  size  is  affected  by  rotating  through  the  inclination  and  declina¬ 
tion.  The  net  effect  of  each  is  influenced  by  the  difference  in  right  ascension,  —  Q), 

“sandwiched”  between  them,  since 
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The  effects  of  each  rotation  in  the  sequence  of  equation  7.6  on  the  apparent  size  and 
shape  of  a  formation  are  graphically  depicted  in  figure  7.5.  The  effect  of  each  rotation 
on  the  intensity  map  is  predictable.  Simple  rotations  in  the  aperture  rotate  the  reference 
intensity  map.  A  decrease  in  aperture  size  along  any  axis  causes  an  expansion  of  the 
reference  intensity  map  along  the  axis  perpendicular  to  the  change,  figure  7.6. 
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Elliptical  Reference  Aperture 
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To  exaggerate  the  effects  of  each  rotation,  consider  a  formation  at  v  =  40°, 
i  =  30°,  and  =  90°  rotated  to  view  a  notional  star  at  <f>  =  30°  and  S  =  —30°. 


Elliptical  Reference  Aperture 


nu  =  40 


nu,  i  =  30 


nu,  i,  (90+RA-AN)  =  30  nu,  i,  AN,  (90+RA-AN),  dec  =  -30 


xi  (rad) 


As  the  elliptical  aperture  is  rotated  and  contracted,  the  effective  intensity 
map  is  directly  affected. 


Performing  the  actual  W*  rotations  for  Andromeda  and  Regulus,  figure  7.7, 


Elliptical  Reference  Aperture  Viewed  from  Andromeda  Viewed  from  Regulus 


Figure  7.7  Aligning  the  iu  axis  of  the  filled  reference  aperture  with  the  look  angle  to 
the  inertial  target  necessarily  warps  the  regular  filled  elliptical  aperture  map. 
Andromeda  is  less  than  50°  off  of  the  centerline  of  the  formation,  and  with 
v  =  40°  the  contraction  is  actually  seen  mostly  in  the  formation’s  major  axis. 
Regulus  is  actually  92°  away  from  the  centerline  (truly  awful  geometry). 


and  the  two  stars  in  Orion,  figure  7.8,  show  the  practical  applications  of  equation  7.18. 


Loci  in  Formation  Frame  Seen  from  Betelgeuse  Seen  from  Rigel 


-50  0  50  -50  0  50  -50  0  50 


u  (wavelengths)  u  (wavelengths)  u  (wavelengths) 

Elliptical  Reference  Aperture  Viewed  from  Betelgeuse  Viewed  from  Rigel 


Figure  7.8  Both  Betelgeuse  and  Rigel  are  less  than  30°  off  of  the  center  axis  of  the  forma¬ 
tion.  Rigel  lies  below  the  Celestial  Equator,  and  shows  a  greater  contraction 
of  the  effective  aperture  size  (wider  peak  in  I0(£,  v) Rigel)- 
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7.5  Optimized  Formations  for  Inertial  Targets 


A  simple  planar  formation  allows  the  optimization  to  select  a  set  of  formation  slopes, 
rnH  and  nH  used  for  all  receivers,  and  both  a  pH  and  dH  for  each  receiver  except  for  the 
one  at  the  origin.  This  results  in  a  problem  with  (Receivers  —  1)  x  2  +  2  degrees  of  freedom. 
Starting  with  a  five  receiver  formation,  four  receivers  must  be  placed  to  provide  a  set  of 
ten  baselines.  Optimizing  such  a  formation  to  view  Andromeda  is  shown  in  figure  7.9. 


Convergence 


Distribution  of  10  Initial  Baselines 


Loci  Viewed  from  Andromeda 
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Ref.  Aperture  viewed  from  Andromeda 


Agg  Map  (pi/2)  4  Receivers 


Figure  7.9  Andromeda  is  only  slightly  off  of  the  formation  centerline,  so  the  reference 
aperture’s  size  is  only  slightly  contracted.  The  initial  v  of  the  formation  has 
the  greatest  effect  on  the  reference  aperture  by  rotating  it. 


Keeping  the  same  formation,  but  optimizing  for  Betelgeuse,  Rigel,  and  Regulus  uses  the 
same  solution  technique,  with  reference  apertures  adjusted  for  each  target,  figures  7.10, 
7.11,  and  7.12. 
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Figure  7.10 


Figure  7.11 


Convergence 
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Ref.  Aperture  viewed  from  Betelgeuse 


Agg  Map  (pi/2)  4  Receivers 


Betelgeuse  also  has  excellent  viewing  geometry  from  the  example  formation, 
showing  only  a  slight  rotation  and  contraction  in  the  reference  aperture. 


Convergence 


u  (wavelengths) 
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Ref.  Aperture  viewed  from  Rigel 


Agg  Map  (pi/2)  4  Receivers 


Rigel  has  fairly  good  viewing  geometry  from  the  formation.  Note  how  the 
optimized  formation  has  been  tilted  with  an  mH  =  9°. 
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Convergence 


Loci  Viewed  from  Regulus 
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Ref.  Aperture  viewed  from  Regulus 
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Figure  7.12  The  poor  viewing  geometry  for  Regulus  is  apparent  both  in  its  reference 
aperture  and  the  selection  of  initial  baseline  orientations.  Note  that  J image 
is  still  quite  small,  indicating  the  formation  is  being  compared  to  an  appro¬ 
priate  reference  aperture. 


Adding  another  receiver  simply  adds  two  more  degrees  of  freedom  to  the  problem. 
Predictably,  the  formations  not  only  mimic  the  filled  apertures  better  with  more  receivers, 
but  the  J image  function  is  decidedly  smaller  for  each  final  formation  (again  suggesting  the 
optimization  is  comparing  appropriate  maps),  figures  7.13  through  7.16. 
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Figure  7.13 


Figure  7.14 
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Ref.  Aperture  viewed  from  Andromeda 


Loci  Viewed  from  Andromeda 
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Optimizing  about  the  same  reference  orbit  with  six  receivers  for  Andromeda 
is  a  12  degree  of  freedom  problem. 


Convergence 
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This  is  the  same  six-receiver  (five  movable)  formation  optimized  for  Betel¬ 
geuse. 
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Figure  7.15 


Figure  7.16 
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This  is  the  same  six-receiver  (five  movable)  formation  optimized  for  Rigel. 


Convergence 
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This  is  the  same  six-receiver  formation  optimized  for  Regulus.  Even  with 
its  atrocious  geometry,  the  optimization  works  well  to  yield  low  values  of 
Jimage- 
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Further  optimization  problems  for  larger  formations  of  7,  9,  and  11  receivers  are 
included  in  appendix  F.  All  show  the  same  trend  as  the  idealized  planar  formation  built 
at  the  end  of  chapter  VI.  Increasing  the  number  of  elements,  N,  in  the  formations  yields 
higher  fidelity  synthesis  of  the  filled  aperture  performance,  but  with  diminishing  returns 
as  N  is  increased. 


5  6  7  8  9  10  11 

Receivers 

Figure  7.17  As  the  formation  is  enlarged,  its  ability  to  synthesize  the  output  of  a  filled 
aperture  increases.  The  aggregate  error  as  compared  to  each’s  filled  aperture 
decreases. 


7. 6  Summary 

This  chapter  expanded  the  models  embedded  in  the  metrics  of  chapter  VI  to  include 
the  distortions  in  both  the  apparent  receiver  loci  and  the  reference  aperture  when  each  is 
viewing  a  target  that  is  not  directly  above  or  below  the  formation  or  aperture.  This  distor- 
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tion  in  loci  was  shown  to  be  applicable  to  the  filled  reference  aperture  yielding  an  analytic 
solution  for  its  intensity  map.  The  technique  yielded  not  only  the  generalized  intensity 
map  for  warped  apertures  but  a  solution  technique  that  may  be  applied  to  most  other 
navigation  or  mapping  conventions  used  to  model  such  systems.  The  distorted  receiver 
loci  and  reference  apertures  in  the  example  problems  continued  to  show  the  expected  im¬ 
provement  in  performance  (decrease  in  mean  aggregate  error)  as  the  number  of  formation 
elements  increased. 
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VIII.  Conclusions  and  Recommendations 


8.1  Summary 

This  investigation  creates  metrics  based  on  models  of  orbiting  distributed  apertures 
performing  interferometry  missions.  The  metrics  allow  the  geometry  of  formations  to  be 
refined  based  on  an  optimization  process.  (Finding  the  true  global  optimum  is  a  more 
extensive  stochastic  process.)  The  metrics  perform  a  comparison  between  the  image  (in¬ 
tensity  map)  generated  by  the  distributed  formation  elements  at  any  time  and  the  image 
theoretically  generated  by  a  filled  aperture  of  size  comparable  to  the  formation. 

Analysis  of  the  nature  of  satellite  formations  reveals  that  for  any  accepted  model, 
a  non-dispersing  formation  limits  design  degrees  of  freedom  to  a  maximum  of  four  per 
receiver.  One  of  the  simplest  (and  therefore  the  lowest  fidelity)  models,  Hill’s  equations, 
is  selected  as  the  primary  model  for  analysis.  The  relatively  short  duration  dwell  times 
for  image  generation  allow  for  this  simplification,  and  the  choice  of  imaging  model  is 
independent  of  the  model  (certainly  of  higher  order)  to  be  used  for  long-term  control  and 
station- keeping.  From  Hill’s  equations,  the  relative  paths  for  each  receiver  element  can  be 
completely  specified  in  terms  of  the  parameters  pH.  0H,  rnH,  and  nH. 

Analysis  of  the  nature  of  interferometers  reveals  that  the  relative  length  and  orienta¬ 
tion  of  the  line  between  each  pair  of  receivers  is  all  important.  Each  baseline  produced  by 
a  receiver  pair  constitutes  a  measurement  of  a  single  spatial  frequency.  (Filled  apertures 
sample  all  spatial  frequency  orientations  and  magnitudes  out  to  the  limit  of  their  diam¬ 
eters.)  Sampling  higher  spatial  frequency  responses  yields  high  image  resolution,  but  at 
the  cost  of  large  blurring  sidelobes.  Lower  spatial  frequency  samples  reduce  the  blurring 
sidelobes  at  the  cost  of  reducing  the  clarity  of  fine  detail.  Selection  of  an  optimal  set  of 
spatial  frequency  magnitudes  and  orientations  is  shown  to  be  a  tradeoff  between  these 
two  effects.  The  actual  generation  of  an  image  requires  sampling  these  selected  spatial 
frequencies  and  performing  their  inverse  Fourier  transform  to  isolate  the  image  intensity 
map,  I (£,77),  in  terms  of  the  spatial  variables  (£,77)  across  the  image  plane. 

Incorporating  time  averaged  errors  into  the  metric  offers  the  promise  of  generating 
solutions  for  which  aggregate  errors  are  near  minimum  for  any  imaging  start  time  at  the 
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cost  of  requiring  consideration  of  the  loci  shapes.  Since  the  formation  elements  do  not  trace 
out  simple  circular  paths,  the  circular  aperture  references  of  previous  works  are  discarded 
in  favor  of  aperture  shapes  that  more  precisely  encompass  all  receiver  loci.  A  new  solution 
for  the  output  of  an  elliptical  aperture  is  derived  as  the  shape  correctly  matching  these 
loci. 

Incorporating  the  distorting  effects  of  viewing  formation  paths  from  off-center  direc¬ 
tions  is  further  investigated.  The  receiver  paths  are  mapped  to  the  inertial  reference  frame 
to  mimic  their  apparent  motion  as  viewed  from  inertial  targets.  This  mapping  depends  on 
the  choice  of  reference  orbit  for  the  formation  and  on  the  inertial  location  of  the  target.  In 
all,  five  angular  variables,  u.  i,  and  f I  from  the  orbit,  and  $  and  5  from  the  inertial  target 
are  incorporated  into  mapping  the  loci.  Beyond  the  individual  loci,  the  filled  reference 
aperture  also  requires  remapping  to  the  inertial  frame.  This  is  accomplished  through  a 
generalization  of  the  derived  elliptical  aperture  to  yield  a  solution  to  a  warped  elliptical 
aperture  that  encompasses  all  possible  receiver  loci  and  no  additional  unreachable  aperture 
areas. 

The  problem  as  posed,  even  for  the  simplest  one-dinrensional  formation,  is  non- 
convex.  For  reliable  solution  approximations  this  necessitated  the  use  of  global  search 
techniques  that  do  not  get  trapped  by  local  minima.  Genetic  algorithms  are  selected  over 
simulated  annealing  or  random  walk  techniques  for  their  robustness  and  degree  of  user  con¬ 
trol.  A  genetic  algorithm  is  constructed  to  generate  generic  populations  of  chromosomes, 
each  yielding  a  dimensionless  coefficient  in  the  range  [—1.0 ...  1.0].  These  coefficients  are 
readily  adapted  to  the  dimensionality  of  any  and  all  of  the  selectable  Hill’s  parameters 
for  each  optimization  problem.  The  search  routines  maintain  a  healthy  “population”  of 
the  most  promising  solutions,  performing  crossovers  and  mutations  between  population 
members  to  generate  new  prospective  members.  The  parameters  defined  by  each  new  set 
of  members  is  subjected  to  a  discrete  time  simulation  that  propagated  each  receiver  out 
in  time,  calculates  the  aggregate  aperture  performance  at  each  time  step,  compares  this 
performance  to  an  appropriately  shaped  reference  aperture,  and  yields  a  time  averaged  ag¬ 
gregate  error.  After  an  empirically  determined  sufficient  number  of  generations,  the  best 
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approximation  of  the  minimum  error  yields  the  set  of  initial  conditions  for  the  formation 
to  give  control-free  performance  with  near-minimal  error. 

The  validity  of  the  constructed  models  is  supported  by  the  behavior  of  J image-  The 
objective  function  approaches  zero  (not  a  fixed  finite  value)  when  the  reference  aperture  is 
shaped  to  precisely  cover  the  receiver  loci  and  the  number  of  receivers  gets  arbitrarily  large. 
This  indicates  that  the  intensity  maps  used  for  I0(£,  77)  or  for  the  inertial  targets  I0(£>7?)* 
are  indeed  correctly  shaped  to  encompass  all  possible  baselines  reachable  by  receivers  and 
none  that  were  not.  Setting  up  the  search  to  synthesize  a  reachable  ideal  has  several 
benefits.  The  genetic  algorithm  has  a  set  of  initial  conditions  accessing  all  baselines  from 
which  to  select  those  paths  sampling  the  most  critical  spatial  frequencies.  It  also  avoids 
expending  computation  effort  finding  conditions  to  synthesize  areas  of  an  aperture  whose 
spatial  frequency  components  are  unreachable  by  the  formation. 

Refining  the  formation  geometry  for  a  solution  with  near-minimum  time-averaged 
error  yields  a  set  of  defining  formation  parameters  (Hill’s  parameters,  here).  These  could 
just  as  easily  have  been  four  parameters  from  any  of  the  models  in  chapter  2.  When 
adjusted  to  these  parameters,  a  formation’s  control-free  dynamics  drove  receiver  loci  to 
yield  the  highest  fidelity  image  for  any  start  time,  rather  than  a  single  critical  point  in 
time.  The  global  search  methods  employed  here  could  be  replaced  with  a  comparable  global 
search  technique,  but  the  fidelity  of  the  reference  map  and  the  utility  of  the  minimum  time- 
averaged  error  remain. 

8.2  Research  Contributions 

This  research  advances  the  science  of  modelling  the  performance  of  satellite  for¬ 
mations  performing  interferometry.  These  advanced  models  facilitate  refining  satellite 
formation  geometries  principally  through  their  use  in  a  performance  metric.  The  search 
technique  employing  this  metric  is  a  simple  genetic  algorithm.  The  highlights  of  these 
advancements  are  as  follows: 
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i.  This  research  formalizes  the  model  for  building  an  intensity  map  generated  by  an 
asymmetric  aggregate  aperture  of  discrete  receivers,  Ip(£,  rj,  t),  with  the  permutation 
operator,  P^,  of  equation  4.12. 

ii.  The  generalized  second-order  metric  of  equation  4.17  includes  model  variations  over 
a  finite  dwell  time.  This  generalized  metric  results  in  formations  whose  control-free 
motion  yields  reduced  imaging  errors  for  any  arbitrary  start  time,  but  requires  a 
more  sophisticated  reference  aperture. 

iii.  Discarding  the  mathematically  simple  circular  reference  aperture  in  favor  of  more 
complex  shapes  includes  all  available  receiver  loci  and  none  that  are  unreachable 
due  to  the  formation  dynamics.  This  requires  a  new  solution  for  elliptical  apertures, 
equation  5.23. 

iv.  The  receiver  loci,  as  viewed  from  an  inertial  target  suffer  some  distortion  that  affects 
the  formation’s  performance.  This  requires  an  enhancement  to  Hill’s  equations  to 
propagate  out  the  receiver  paths  as  viewed  from  an  inertial  location,  equation  7.6. 

v.  The  analytic  solution  for  an  elliptical  aperture  is  generalized  to  model  the  distorting 
effects  of  non-ideal  viewing  geometries  in  terms  of  the  Keplerian  elements  of  the 
reference  orbit  and  the  inertial  locations  of  targets,  equation  7.18.  This  warps  the 
known  elliptical  aperture  solution  to  encompass  all  possible  receiver  paths  as  viewed 
from  each  inertial  target. 

vi.  The  method  used  to  generalize  the  filled  aperture  solution  is  a  very  generic  technique 
for  mapping  aperture  distortions  through  the  PSF  to  the  image  plane.  It  may  be 
adapted  to  model  the  distortions  of  viewing  geometries  described  by  essentially  any 
mapping  convention. 

vii.  The  genetic  algorithms  used  were  are  not  advanced,  but  do  illustrate  the  problem 
formulation  for  selecting  the  defining  parameters  of  simplified  formations.  This  for¬ 
mulation  and  algorithm  application  can  be  adapted  to  more  advanced  algorithms  or 
to  higher  fidelity  models,  such  as  the  four  Delaunay  elements  of  Schaub  and  Alfriend 
or  the  four  modal  variables  of  Wiesel,  with  minimal  modification. 
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8.3  Recommendations  for  Future  Research 

As  posed,  the  optimization  problem  is  notoriously  non-convex,  necessitating  the  use 
of  computationally  intensive  global  search  techniques.  The  architecture  and  tools  provided 
by  Optimal  Synthesis®  in  their  Genetic  Search  Toolbox®  are  user-friendly  (after  a  fashion) 
but  carry  the  heavy  computation  load  of  the  MatLab  environment  (acceptable  only  on 
very  fast  workstations).  Though  larger  formations  should  be  expected  to  require  longer 
computation  times,  the  durations  become  impractical  for  formations  with  more  than  a 
dozen  receivers.  More  efficient  genetic  algorithms  exist  and  should  be  readily  adaptable 
to  the  problem.  With  more  efficient  computation,  further  assumptions  on  the  formation 
geometry  may  be  relaxed,  principally  the  planar  assumption.  By  allowing  the  full  four 
degree  of  freedom  placement  of  each  receiver’s  relative  orbit,  the  true  optimum  may  be 
more  closely  approximated. 

The  formation  model  used  throughout  this  work,  Hill’s  model,  was  the  simplest  and 
therefore  the  lowest  fidelity  model  available.  The  Wiesel-Floquet  solution  (Appendix  C) 
offers  modelling  accuracies  several  orders  of  magnitude  better  than  Hill’s.  Including  such 
sophistication  into  the  optimization  model  is  intriguing,  but  with  the  understanding  that 
this  level  of  fidelity  would  really  only  be  necessary  for  systems  with  very  long  imaging  times 
or  systems  using  very  short  wavelengths  and  requiring  the  accompanying  high-precision 
receiver  positioning. 

The  non-convexity  of  the  problem  may  well  be  removed  through  the  application  of 
new  convex  hull  techniques.  A  significant  amount  of  research  effort  was  performed  (but  not 
included,  though  the  preliminary  modelling  was  promising)  to  formulate  the  optimization 
problem  with  different  logic  than  that  used  here.  Rather  than  optimize  each  formation  in 
a  growing  series,  the  problem  may  be  approached  by  weighting  the  contributions  of  a  large 
number  of  receivers.  Minimizing  the  sum  of  the  weights  (an  L±  metric)  while  constraining 
the  aggregate  aperture  performance  to  an  arbitrary  tolerance  holds  the  promise  of  zero 
or  near-zero  weights  on  many  of  the  receivers.  Thus  the  solution  becomes  a  matter  of 
reducing  a  formation  size  by  eliminating  the  least  productive  receivers. 
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Appendix  A.  The  Clohessy-  Wiltshire  Solution 
The  Clohessy-Wiltshire  approach  to  describing  the  dynamics  of  clustered  satellites  in  a 


local  coordinate  frame  began  with  choosing  a  circular  orbit  of  radius  Ra  as  the  local  origin. 


with  mean  motion  uj  =  y  Though  Clohessy  and  Wiltshire  originally  chose  the  ey 
axis  in  the  radial  direction  and  the  ex  axis  in  the  negative  velocity  direction  [5],  the  more 
conventional  approach  sets  ex  in  the  radial,  ey  in  the  velocity  direction,  and  ez  normal  to 
the  orbit  plane  [49] .  (Forgive  the  volume  of  algebra  presented  hereafter,  but  it  is  worthwhile 
to  see  the  extensive  approximations  necessary  to  yield  a  linear  system.) 


Figure  A.l  The  Clohessy  and  Wiltshire  problem  is  more  conventionally  addressed  with 


ex  in  the  radial  direction  and  ey  in  the  velocity  direction. 


In  this  reference  frame,  the  radius  vector  to  a  cluster  satellite  at  some  location  (x,  y,  z )  in 
the  local  frame  is 


r  =  (x  +  Ro)e  x  +  yey  +  zez. 


(A.l) 


Since  the  local  reference  frame  rotates  with  u,  the  velocity  vector  is 


(A.2) 


and  the  acceleration  vector  is 


A-l 


r=  [(x  -  uy)ex  +  (y  +  ux)ey  +  zez]  +  [u(x  -  toy)ey  -  uj(y  +  uj(x  +  R0))ex\  ■  (A-3) 

f=(x  —  2uiy  -  uj2(x  +  R0))ex  +  (y  +  2w.x  -  ui2y)ey  +  zez 

In  the  absence  of  control  forces,  gravity  is  the  only  acceleration.  Therefore 


-n((x+Ro)ex+yey+zez) 

((i+B„)2+!/2+z2)2 


(A.4) 


The  denominator  here  is  heinously  non-linear,  but  it  may  be  expanded  using  the  binomial 
theorem. 


((x  +  Ro)2  +  y2  +  z2)  2  =  (R20  +  2 R20x  +  x2  +  y2  +  z2)  2 

((*  +  +  j,2  +  ,2)f  =  .R’  [l  +  (g  +  ff  +  ^  +  |)] 1  (A.5) 

((*  +  Rof  +  f  +  *?)§  =  JJ3  [1  +  |  (g.  +  g  +  ||  +  g)  +  H.O.T. 

Assuming  small  components  in  the  local  coordinate  system,  |^,  and  may  be  as¬ 
sumed  «  0,  leaving 

q 

{{x  +  R0)2  +  y2  +  z2)~2  *tRl  1  +  —  =R30  +  3R20x  (A.6) 

rC0  _ 

which  is  closer  to  linearity.  Inserting  this  into  the  r  expression, 

-  -y((x  +  R0)ex  +  yey  +  zez) 

r“ - iiFTvfP - '  <A'7) 

This  expression  for  the  acceleration  due  to  gravity  can  be  equated  component  by 
component  with  equation  A. 3  to  yield  the  equations  of  motion. 
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ex  ■  R3+3R2°x  =  X-  2ujy  -  u2{x  +  Ro) 

:  Ki+zkix  =  y  + 2uj±  -  u2y  (A-8) 

~  —fiz  _  •• 

:  Ki+ZKix  -  Z 

Noting  that  n  =  to2 R°,  a  quick  bit  of  algebra  yields 


x  =  2 uy  +  u?[x  +  R0)  -  ^  f°R0t^ 

^oK  Rq  ) 

y  =  -2“X  +  u2y  -  ^k5-L} 

••  _  u)2R^z 

-  - 


(A.9) 


x  =  2 ujy  +  lu2(x  +  R0 )  -  lo2(x  +  R0)  (r^+[ 
y  =  -2ux  +  uj2y-uj2y[^^j 


Z  =  —UJZ 


Ro 


Ro+3x 


(A.10) 


x  =  2tvy  +  u2(x  +  R0)  (j^) 

y  =  -2^ +  ^1/(7^) 


Z  =  —UJZ 


Ro 


R0+3x 


(A- 11) 


For  small  displacements  in  the  local  coordinate  frame,  it  is  possible  to  linearize  further. 
Assuming  x  +  R0  ~  RQ,  and  3x  +  R0  «  Ra, 


0 

x  =  2 ujy  +  3uj  x 

(A. 12a) 

y  =  —2  lux 

(A.  12b) 

z  =  —oj2z 

(A. 12c) 

which  are  the  classical  Clohessy-Wiltshire  equations  of  motion. 

By  employing  normalized,  canonical  units,  (u  =  1,  /r®  =  1)  it  is  possible  to  array  the 
state  vector  to  set  up  a  linear  system.  Adding  in  thrusts  along  each  axis  as  Tx,  Ty,  and  Tz, 
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(A.13) 


This  has  the  familiar  form  of  a  linear,  constant-coefficient  dynamic  system  [24]. 


x  =  Fx  +  Bu 


(A.  14) 


A-4 


Appendix  B.  Hill’s  Equations 


As  derived  in  Appendix  A,  the  linearized  equations  of  motion  for  a  satellite  formation  in 
a  localized  coordinate  system  are 


X  = 

2  uiy  +  3  uj2x 

(B.la) 

y  = 

—2  ujx 

(B.lb) 

z  = 

—0J2Z. 

(B.lc) 

As  presented  by  Yeh  and  Sparks  [56],  these  equations  can  be  parameterized  in  terms 
of  six  arbitrary  constants:  a,  b,  p,  m,  n,  and  9.  Summarizing  the  derivation  by  Irvin  [15], 
these  equations  can  be  solved  in  terms  of  these  constants  by  first  performing  the  Laplace 
transformation 


s2Y(s)  —  sxQ  —  x0  =  2w[sY(s)  —  y0\  +  3  lo2X(s) 

s2Y {s)  -  syQ  ~y0=  -2 w[sl(s)  -  x0]  (B.2) 

s2Z(s)  —  sz0  —  z0  =  lu2Z(s), 

where  the  subscripted  “0”  values  are  the  components  of  initial  position  and  velocity.  Col¬ 
lecting  these  initial  conditions,  these  equations  may  be  expressed  in  matrix  form  as 


s2  —  3u>2  —2  ws 

0 

'  X(s)  ' 

sx0  +  Xo  -  2 ivy0 

2  cos  s 2 

0 

Y(s) 

= 

sy0  +  Vo-  2c ox0 

1 

o 

o 

s2  +  UJ2 

1 

i _ 

sz0  +  z0 

Solving  for  X(s),  Y(s),  and  Z{s) 
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'  X(s)  ' 
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2  LJ 
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sx0  +  x0-  2u )ya 

s2+lj2 

s(s2+  a;2) 

Y(s) 

_ 

— 2cj 

s2—3uj2 

0 

sya  +  y0-  2wx0 

s(s2+lo2) 

ss(s'2+lj2) 

Z(s) 

0 

0 

1 

SZo  +  Zo 

S2+UJ2  _ 

(B.4) 
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"  X(s)  ' 

sXo+Xo- 2ujy0  ,  2uj(sy0+y0)+2ux0 
s2+u>2  1  s(s2+u>2) 

Y(s) 

_ 

-2oj(sx0+x0-2ujy0)  ,  (s2-3u2)(sy0+y0+2ujx0) 

s(s2+oj2)  s2(s2+oj2) 

Z(s )  _ 

SZo  +  Zo 

S2+(jJ2 

Performing  the  partial  fraction  expansion, 


'  X(s)  ' 
Y(s) 

sXo+Xo—^uyo  |  2uj2y0— 2sy0— 4so;rr0  i  2y0-\-^x0  1 

S2+UJ2  Oj{s2+OJ2)  OJ  s 

4sy0+4y0+8ux0  ,  -2cj2x0+2sx0-4sy0u>  ,  -3y0-6ujx0  ,  u)y0-2x0  1 

s2+u>2  ui(s2+oj'2)  s'2  UJ  s 

Z{s) 

sz0+z0 

s'2+UJ2 

Set  up  for  the  inverse  Laplace  transform  by  collecting  terms, 


'  X(s)  ' 

r  (x„-2-^±^)s  (£)W  2*0+4W*ol  1 

S2+tt)2  S2+W2  OJ  s 

Y(s) 

2(uj)s  |  2(  u,  X°)U}  3u  2y0+4u>x0  1  .  toy0-2x0  1 

s2+u;2  s2+w2  2  o>  s2  '  lj  s 

Z(s) 

(ie-)cj 

SZ0  l  \  u>  / 

S2+0)2  S2+OJ2 

and  perform  the  inverse  transform  [36] 


x{t) 

y(t) 


z(t) 


— 3x0  —  cos(c ot  +  6)  +  —  sin(c ot  +  9)  H — —  +  4x0 

to  )  to  to 


LO 


2  y 

cos (tot  +  6)  +  2  (  — -  +  3x0  )  sin(u;t  +  0) 


3 wt  (  2y0 

-^Ur+4a:”l  + 


to 

uy0  -  2 x0 

LO 


z0  cos(  tot  +  6)  H — -  sin  (tot  +  6) . 

LO 


(B.8a) 

(B.8b) 

(B.8c) 


Note  that  6  is  just  an  initial  phase  angle.  Since  the  initial  conditions  are  not  time- 
dependent,  and  a;  is  a  constant  for  the  system,  it  is  possible  to  define  the  following  con¬ 
stants, 


B-2 


aH 


Pi 


mH 


nH 


tan  0H 


2  yQ  +  4  u)x0 

UJ 

ivy0  -  2x0 
u 

/  \2  ,  f  X°\ 

\Xq  —  ci)  +1  —  ] 

\  UJ  J 

ZqXq  -  z0u?[a  -  Xp) 
+  co2(a  —  x0)2 
z0x0lu  +  z0uj(a  —  x0) 
2(x2  +  co2(a  -  x0)2) 
u>(x0  —  a) 
x0 


(B.9a) 

(B.9b) 

(B.9c) 

(B.9d) 

(B.9e) 

(B.9f) 


These  constants  describe  the  size,  shape,  and  phase  of  the  formation  elements’  relative 
orbits  in  the  local  coordinate  system.  Their  physical  interpretation  within  the  x  —  y  plane 
is  shown  in  figure  B.  The  parameter  rnH  is  the  tangent  of  the  angle  between  the  minor 
axis  and  ex,  and  nH  is  the  tangent  of  the  angle  between  the  major  axis  and  ey. 


Figure  B.l  The  constants  defined  to  parameterize  the  C-W  solution  define  the  relative 
orbit  size  ( pH ),  location  ( aH  and  bH),  and  phase  ( 0H ). 

Using  these  parameters,  Clohessy  and  Wiltshire’s  solutions  may  be  written  as 
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x{t) 

=  pHsin(u>t  +  9h)  +  aH 

(B.lOa) 

y(t) 

,  .  3  LO 

=  2pHcos(ujtJrf)H) - —aHt  +  bH 

(B.lOb) 

z(t) 

=  mHpH  sin (ivt  +  0H)  +  nH2pH  cos {uit  +  0H). 

(B.lOc) 

Obviously,  the  secular  term  in  y(t)  causes  a  problem  for  maintaining  a  closed  path  in 
the  local  coordinate  system.  Therefore,  keeping  a  formation  together  requires  aH  =  0, 
or  setting  the  center  of  the  orbit  on  the  ey  axis.  This  yields  the  initial  constraint  of 
y0  =  —2 lox0.  From  Yeh  and  Sparks  [56],  the  parameterization  yields  the  identities 

(x—aH)2  ,  {y+^faHt-bH)  _  , 

ph  (B.n) 

z(t)  =  mH(x  -  aH)  +  nH{y  +  ^ aHt  -  bH). 

For  a  non-dispersing  formation,  these  become 

x 2  ,  (■ y-bH )2  _  -I 

^  -  (B.12) 
z(t)  =  mHx  +  nH(y  -  bH), 

and  bH  becomes  simply  an  offset  of  the  relative  orbit  along  the  velocity  direction.  Since 
Hill’s  equations  have  a  limited  range  of  validity,  there  is  no  immediate  reason  to  offset  the 
origin  from  the  ellipse  center,  so  formations  also  set  bH  =  0.  Thus  the  degrees  of  freedom 
for  the  orbit  of  each  allowable  element  of  a  formation  have  been  reduced  from  six  to  four, 
in  agreement  with  the  initial  formulations  of  Schaub  and  Alfriend  [37] . 
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Appendix  C.  The  Wiesel-Floquet  Solution 

As  presented  by  Dr.  William  E.  Wiesel  [54],  the  linear,  time-invariant  solution  of  Clohessy- 
Wilshire  can  be  replaced  by  a  linear,  time-periodic  solution  using  the  techniques  of  Floquet. 

C.l  The  Hamiltonian  in  the  C-W  Reference  Frame 

Begin  by  defining  X,  the  three  inertial  coordinates  and  three  inertial  momenta  (per 
unit  mass,  essentially  the  velocity  components)  of  a  formation  element,  as  X7  =  [ X ,  Y,  Z,  P\,  Py,  Pz]  ■ 
Using  normalized,  canonical  units  (similar  to  the  end  of  Appendix  A,  where  Rq  =  1, 

M®  =  1,  and  G  =  1),  all  of  the  zonal  harmonics  may  be  included  in  the  Hamiltonian,  Tt, 
to  yield 

Z2 

—  30— 2 — b  3 
"(C.l) 

where  r  =  \/ X'2  +  Y2  +  Z2.  Additional  conservative  perturbations,  besides  the  zonal  har¬ 
monics  shown,  may  be  added  to  this  Hamiltonian  by  simply  adding  in  their  potential 
functions.  Non-conservative  perturbations  (drag,  etc.)  are  slightly  more  difficult,  but  may 
be  added  by  first  expanding  them  about  the  periodic  orbit  and  then  adding  them  in  terms 
of  the  momentum  states. 

In  inertial  space,  this  orbit  will  not  close  on  itself,  since  it  undergoes  some  nodal 
regression,  H.  Setting  up  a  rotation  matrix  for  a  total  rotation  of  +  Qt  about  the 
inertial  Z  axis  is  simple. 

cos(H0  +  tit)  sin(H0  +  tit)  0 

Hz  =  —  sin(H0  +  Hf)  cos(H0  +  tit)  0  (C.2) 

0  0  1 

This  allows  the  components  of  X  to  be  mapped  into  the  rotating  nodal  reference  frame,  J\f, 
via 


n  -  o  iPx+PY  +  Pz}-Z+Yf3  1  3XX  “  1  \  5XX  -3X  f  +  s^5  1  35X4 
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P' 
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P' 

rz  J 
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.  pz  . 

(C.3) 


With  this  rotation  into  the  nodal  frame,  the  Hamiltonian  grows  an  extra  term.  Using 
V (Z',  r )  to  incorporate  the  zonal  harmonic  terms  from  equation  C.l  rotated  into  the  nodal 
frame,  the  new  Hamiltonian,  TC ,  becomes 


H'  =  \  {p'x  +  PY  +  P'l)  +  { P'xY '  -  P'YX '}  +  V{Z',r).  (C.4) 

In  this  frame,  the  orbit  now  closes  on  itself  and  is  therefore  periodic,  and  has  a  periodic 
solution.  It  is  possible  to  shift  the  origin  out  to  the  local  coordinate  origin  of  the  familiar 
C-W  reference  frame.  This  is  done  by  subtracting  off  the  reference  orbit  state  vector 
M0  =  (X'0,  Y'0,  Z'Q,  P'Xoi  Py0,  P'zo )  in  the  nodal  frame  to  define  y  as  1 


y(t)  = 


Xn 

'  X'  ' 

’  K  ' 

Vn 

Y' 

Y' 
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Z' 

Z'o 

Xn 
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PX  o 

Vn 

P' 

rY 

pi 

r Yo 

Zn 

pi 

L  rz  J 

1 

/t) 

o 

i _ 

(C.5) 


This  is  a  canonical  transformation,  and  reduces  the  Hamiltonian,  7 i." ,  function  to  its 
quadratic  and  higher  order  terms  evaluated  on  the  periodic  orbit  (0  subscript) 


1The  notation  for  coordinates  hereafter  is  a  significant  departure  from  that  used  by  Dr.  Wiesel  [54],  but 
is  adopted  for  consistency  with  that  used  in  the  C-W  solution  and  Hill’s  equations,  where  the  notation 
(x,  y,  z,  x,  y,  z)  is  reserved  for  coordinates  specifically  in  the  C-W  reference  frame. 
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n ill 
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343^ 


(C.6) 


where  the  repeated  Greek  indices  simplify  the  notation  for  tensor  multiplication,  and  are 
summed  from  one  to  six.  The  and  terms  are  from  the  fully  symmetric  second, 

third,  and  higher  order  tensors  of  partial  derivatives  with  respect  to  the  nodal  vector,  J\f. 

One  final  rotation  is  necessary  to  align  these  axes  to  the  familiar  radial,  velocity,  and 
orbit  normal  axes  of  C-W.  Using  f  and  v  for  the  nodal  frame  position  and  velocity  vectors, 
define  the  rotation  matrix  1Z  by 


A 

r 

Ul 

M 

/\  /\ 

U2 

Us  x  Ul 

r  X  v 

U3 

|rx  u| 

(C.7) 


and  finally  write  the  new  position  and  velocity  coordinates  in  the  standard  C-W  frame, 
T’(t),  with 


X 

0  0  0 

y 

0 

0 

0 

& 

Vn 

z 

0  0  0 

Zn 

X 

0  0  0 

y 

000  n 

Vn 
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000 

Zn 

(C.8) 


This  is  also  a  canonical  transformation,  with  a  generating  function 


X 

X 
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■  n 

yn 

= 

y 
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Zn 

z 

z 

(C.9) 


with  which  the  Hamiltonian  can  be  transformed  term  by  term  into  the  X  variables  to  yield 
1C.  The  quadratic  term,  1C 2,  looks  like  this 
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K2  =  hKByc,ys+d-& 


X 

X 

1C2  =  hK^lKpsX'yXs  + 

y 

■  nnT 

y 

z 

z 

(C.10) 


where  TZ^  is  the  block  diagonal  rotation  matrix  from  equation  C.8.  Higher  order  terms 
are  even  simpler;  the  third  order  term  is 


£3 


n[2l  n xaxpx7  • 


g  I  ^  A  err  /VA  a  /vcr/3  ,VT7 


(C.ll) 


C.2  The  Floquet  Solution 

Though  the  Hamiltonian  derived  for  the  C-W  reference  frame  is  truly  an  infinite 
sum,  approximating  it  by  only  the  second  order  term,  1C 2,  not  only  preserves  inclusion  of 
all  the  zonal  harmonics,  but  yields  a  matrix  linear  system.  Using  Hamilton’s  equations  of 
motion  [11]  on  !C(x,y,  z,x,y,  z),  it  is  possible  to  generate  the  plant  matrix  from 


™  _  _8K  •;  _  _  <WC  -  _  _dK. 

dx  ’  &  dy'^  dz 

such  that  the  equations  of  variation  may  be  written 


dX  =  A  {t)dX 
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(C.13) 


Their  solution  [8],  the  state  transition  matrix  <b(f, 0),  may  also  be  propagated  via  the 
equations  of  variation 

X(t)  =  $(t,0)A(0)  (C.14) 
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$(t,0)  =  A(t)$(*,0). 


(C.15) 


Since  the  A(t)  matrix  is  periodic,  A (t  +  r)  =  A (t),  the  state  transition  matrix  for  the 
reference  orbit  (whether  there  is  a  satellite  at  the  origin  or  not)  may  be  integrated  over 
one  orbit  period,  r,  which  yields  the  monodromy  matrix,  $(r,  0).  This  state  transition 
matrix  may  be  factored  using  the  techniques  of  Floquet  [25]  as 

$(t,  0)  =  F(t)eJfF(0)-1,  (C.16) 

where  F(t)is  a  periodic  modal  matrix,  and  Jis  a  Jordan  normal  form  matrix  of  system 
frequencies  described  by  Poincare  exponents,  of  the  form 


0  0  0  0  0 

-wi  0  0  0  0  0 

0  0  0  1  0  0 

0  0  0  0  0  0 

0  0  0  0  0  1 

0  0  0  0  0  0 


(C.17) 


where  the  one  non- degenerate  frequency  uj\  is  the  frequency  of  the  advance  of  perigee. 
The  frequencies  may  theoretically  be  determined  from  the  eigenvalues  of  the  monodromy 
matrix,  \  via 


\  =  e^r 

(C.18a) 

In  A  i 

LOi  — 

T 

(C.18b) 

Using  this  relation,  and  substituting  equation  C.16  into  equation  C.13  yields 

F(t)  =  A(t)F(f)  -  F(f)J,  (C.19) 

from  which  F (t)  may  be  found.  However,  this  is  a  difficult  integration,  since  there  are  two 
integrals  of  the  motion,  two  pairs  of  the  Poincare  exponents  are  zeros.  Therefore  all  but  two 
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of  the  eigenvectors  of  F(f)  are  generalized  eigenvectors.  Wiesel  devised  a  mildly  convoluted 
method  for  directly  determining  these  eigenvectors  (which  most  software  packages  utterly 
fail  to  accomplish). 


Label  the  first  two  eigenvectors  fi  and  £2,  corresponding  to  the  non-zero  eigenvalues. 
Two  more  eigenvectors  may  be  directly  calculated.  Use  k  as  the  unit  vector  in  the  inertial 
Z  direction,  and 


f3 


-p  dNo 
dt 


’  X'  ’ 

7 Zk  x 

Y' 

7 Zk  x 
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z 

(C.20) 


The  remaining  eigenvectors,  f4  and  fg,  must  satisfy 


-  I)  U  =  rf3 
($-I)f6  =  rf5, 

but  it  is  not  possible  to  directly  calculate  f4  and  fe  with  (<h  —  I)-1  since  (<h  —  I)  has  four 
zero  eigenvalues.  Noting  that  the  generalized  eigenvectors  are  undetermined  to  within  the 
addition  of  a  multiple  (7)  of  one  of  the  normal  eigenvectors,  the  remaining  eigenvectors 
must  satisfy 

=  f4  +  O'fa 

fe  =  f6  +  7f5- 

Finally  the  orthogonality  of  the  eigenvectors  may  be  used  to  resolve  the  ambiguity  since 


(C.22) 


f4  •  f3  =  fe  •  f3  =  0 

f4  •  f5  =  f6  •  f5  =  0. 


(C.23) 
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With  the  F (t)  matrix  calculated,  the  position  vector  of  a  satellite  may  be  determined 
in  the  C-W  reference  frame  by  inserting  equation  C.16  back  into  C.14 

X(t)  =  F(t)eJtF(0)_1*(0).  (C.24) 

The  position  vector  in  inertial  coordinates  may  be  determined  by  rotating  the  C-W  refer¬ 
ence  frame  vector  to  the  nodal  coordinate  frame  using  1Z(-2',T,  adding  in  the  reference  orbit, 

(o\rr 

N0,  and  rotating  the  sum  into  the  inertial  frame  with  7 Zyz 

I{t)  =  1Z(z)T  ( M0{t )  +  ft(2)TF(i)eJtF(0)_1*(0))  (C.25) 

C.3  The  Modal  Variables  Solution 

The  raw  Floquet  solution  of  the  last  section  makes  physical  interpretation  difficult. 
Therefore,  Dr.  Wiesel  proposed  the  further  change  of  variables  from  the  variational  equa¬ 
tions  of  equation  C.13  to  the  modal  variables  Z.  Define 


y(t)  =  TL(2\t)F(t)Z  =  E  (t)Z,  (C.26) 

where  Z  is  a  six-vector  of  the  Floquet  modal  amplitudes,  Z 1  =  (z\,  Z2,  Z3,  Z4,  Z5,  zq),  and 
E(t)  is  the  Floquet  modal  matrix  rotated  to  the  nodal  frame.  Substituting  this  into 
Hamilton’s  equations  of  motion,  written  for  hi! ,  yields 


ftz  =  (e_1ae  -  e_1e)  z 
ftz  = 


(C.27) 


where  J  is  the  Jordan  normal  form  matrix  from  equation  C.17.  Dr.  Wiesel  showed  [52] 
[53]  that  with  J  arranged  in  this  manner,  Z4  and  z%  must  be  set  to  0  to  prevent  formation 
dispersion,  again  in  agreement  with  the  four  degree  of  freedom  conclusion  of  Schaub  and 
Alfriend  [37]. 
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Appendix  D.  Spatial  Frequency  Measurements  in  Interferometery 

This  derivation  summarizes  Thompson’s  [46]  description  of  how  coordinates  in  the  {u  v  w) 
reference  frame  are  both  unitless  physical  lengths  and  spatial  frequency  variables  (as  could 
be  inferred  from  their  appearance  in  V(u,v)). 

The  basic  operation  of  an  interferometer  analyzes  an  incoming  wavefront  from  a 
distant  source  as  received  by  two  receivers  with  separation  D.  Take  the  geometry  from 
figure  3.4,  and  assume  for  the  moment  that  the  antennas  have  perfect  beam  patterns,  so 

F(D)  =  R(D). 


Figure  D.l  Incident  wavefronts  from  a  target  source  at  6  arrive  at  the  two  receivers  at 
times  that  differ  by  r  =  Ds™° . 


The  output  of  the  correlator  at  any  time,  t,  should  be  a  simple  function  of  the  center 
frequency  is0,  the  baseline,  and  the  angle  to  the  source. 


R(D,0,  t)  =  2  sin[27rzz0t]  sin 


27 TU„  t  — 


=  COS  27TZA, 


D  sin  6 


D  sin  6 


—  cos(47TZ70t)  cos  (  2nu, 


D  sin  i 


—  sin(47rt'0f)  sin  2ttu, 


D  sin  6 


A  correlator  also  serves  to  filter  out  noise.  Since  the  vQ  is  far  greater  than  the  frequency 
of  6  (the  angular  rate  of  physically  moving  the  receivers),  the  high  frequency  portions  of 
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R(D,0,  t)  are  easily  filtered  out,  leaving1 


R(D,  9)  ~  cos  (  - - —  sm  ( 


/27rT>  . 

=  cos  — —  sm  9 


An 


(D.2) 


where  Ac  is  the  wavelength  of  the  center  frequency.  Since  a  correlator  both  multiplies  and 
time-averages  the  received  signals  (voltages,  V(t)),  the  response  over  a  finite  time  T  is 


R(D,  9)  =  R(r)  =  J  V(t)V  ( t  -  r)  dt. 


(D.3) 


In  the  limit, 

R(r)  =  lirn  [  V{t)V{t-r)dt  (D.4) 

1  — >oo  ZjI  J —T 

which  Thompson  terms  an  autocorrelation  function  (or  is  directly  proportional  to  Rohlfs  [35] 
“mutual  coherence”  function,  equation  3.6). 

Rather  than  assume  a  constant  field  strength  from  the  source,  as  in  the  derivation  of 
section  3.2.1,  allow  it  to  vary  by  frequency.  At  any  frequency,  u,  the  field  strength,  U(i')  is 
directly  related  to  the  response  function  by  the  Wiener-Khinchin  relation  (see  Thompson 
[46]  and  Maybeck  [24]).  The  power  spectral  density,  \U(v)\2,  of  a  signal  is  the  Fourier 
transform  of  its  autocorrelation  function 


\U(u)\2  =  f°°  R{T)e~2i™TdT 

'  V  '  J-°°  V  '  (D.5) 

R(t)  =  \U{v)\2e2i™Tdv 

Therefore  the  interferometer’s  output  as  a  function  of  9  (or  r)  is  the  Fourier  transform  of 
the  power  spectral  density  of  the  signal.  This  yields  the  general  result  that  the  response 
function  R(r)  shows  fringes  whose  spatial  frequency  period  varies  inversely  with  at  any 
signal  frequency  (just  as  in  equation  3.15). 

In  practice,  to  use  an  interferometer  to  map  a  target,  the  delay  time  r  is  calculated 
for  a  nominal  direction  90  at  the  source  center.  (This  may  also  be  tuned  to  a  desired 
shorter  delay  by  adding  in  an  instrument  delay.)  For  some  small  angular  displacement  9' 

1  Rohlfs  [35]  shows  that  this  approximation  is  valid  for  correlators  designed  to  either  multiply  or  add  the 
signals. 
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from  0o,  the  correlator  output,  equation  D.2,  for  0  =  0o  —  O'  becomes 


R(D,0)  =  R(r) 


cos 

cos 

cos 


D  sin  0  \ 

A0  ) 


27t—  (sin  0O  cos  O'  —  cos  0o  sin  O' ) 
27T— (sin  0o  —  O'  cos  0o) 

^ O 


(D.6a) 

(D.6b) 

(D.6c) 


It  is  therefore  convenient  to  define  a  value  u  =  D  c^s  60  to  be  the  component  of  antenna 
spacing  orthogonal  to  the  direction  0o  and  £  =  sin#'  ~  O'.  The  (one  dimensional)  output 
fringe  pattern  for  a  coordinate  system  oriented  so  0o  =  0  with  fixed  D  becomes 


R(£)  =  cos(27nx£) 


(D.7) 


The  quantity  u  is  interpreted  as  a  unitless  length  describing  a  component  of  antenna 
spacing,  measured  in  number  of  wavelengths  of  Ac.  It  is  also  a  spatial  frequency  with 
dimensionless  units  of  ^ Jdian  >  since  it  describes  the  transformed  response  of  a  source  at  f 
radians  away  from  0o.  (See  equation  3.18  to  quickly  see  why  u,  £  and  their  two  dimensional 
counterparts  v  and  r]  must  be  dimensionless.) 
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Appendix  E.  A  Genetic  Algorithm 

The  purpose  of  this  research  is  not  to  advance  the  science  of  genetic  algorithms.  However, 
the  heavy  reliance  on  them  to  solve  the  inherently  non-convex  problems  warrants  a  more 
extensive  discussion  of  their  implementation. 

The  numeric  portions  of  this  research  were  conducted  in  the  MATLAB  5.0  and  6.0 
environments.  Genetic  search  routines  are  notoriously  computationally  intensive  (solution 
time  ranged  for  24  hours  for  small  formations  to  five  or  six  days  for  the  largest  formations), 
but  the  speed  of  available  workstations  made  the  MATLAB  computational  overhead  an 
acceptable  tradeoff  with  user-friendliness  and  available  genetic  operation  support.  Most 
notably,  the  Genetic  Search  Toolbox  1.0.3  by  Optimal  Synthesis  Inc  provided  a  significant 
set  of  genetic  operation  functions  and  guidance  for  building  stable  search  routines.  Rather 
than  include  a  complete  listing  of  raw  code,  this  appendix  summarizes  the  coding  logic 
and  sequence  of  commands  used  to  build  and  execute  the  genetic  search  “engine.” 


E.l  Defining  the  Initial  Populations 

The  Genetic  Search  Toolbox®  is  designed  to  operate  on  genes  and  chromosomes  as 
character  strings.  Its  basic  set  of  genes  are  defined  as  a  set  of  alpha-numerics,  each  with 
a  numeric  equivalent.  These  may  then  be  converted  to  numerical  values  by  evaluating  the 
numerical  equivalent  of  each  name  with  the  associated  arithmetic  operations.  As  cited  in 
chapter  VI,  the  set  of  single  precision  genes  used  certainly  span  the  solution  space  (but  is 
none  too  clever).  The  actual  coding  looks  more  like  this: 


ul=l 

u2=2 

u3=3 

u4=4 

u5=5 

u6=6 

u7=7 

u8=8 

u9=9 

ul0=0. 

.1 

ul 1=0 . 

.01 

ul2=0 , 

.001 

ul3=0. 

.0001 

ul4=0. 

.00001 

ul5=0. 

.000001 

Using  these  character  string  genes,  it  is  necessary  to  define  an  initial  set  of  chromosomes 
for  each  population.  For  coding  efficiency,  the  same  set  of  chromosomes  is  built  into 
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each  population.  These  then  build  their  own  “diversity”  in  succeeding  generations.  The 
population  building  process  uses  the  gs  new(.  .  .)  command,  in  a  larger  form  of  the 
following  loop. 


for  rr  =  l:NumPops 

mem  id  =  gs  new( [ ’Pop’ ,num2str (rr)] ,  Ju6*ul0+ul5+ul4+ul3+ul2+ull ’ ) ; 

mem  id  =  gs  new( [ JPop’ ,num2str (rr)] ,  Jul5*ul5+u5*ul4+u6*ul3+u7*ul2+u8*ull ’ ) ; 

mem  id  =  gs  new( [’Pop* ,num2str (rr)] ,  5ul5*ul3+ul4/ul2+ul3*ull+ul2/ul0+ull*ul 1 ) ; 

mem  id  =  gs  new( [’Pop’ ,num2str (rr)] ,  5ul5/ul3+ul4*ul2+ul3/ull+ul2/ul0+ull*ul ’ ) ; 

mem  id  =  gs  new( [JPop} ,num2str (rr)] , ’ -ulS’fr^/u^-u^’i'uSA^+ulS’i'M-u^+uS-ulO+u^ ) ; 

mem  id  =  gs  new( [ JPop’ ,num2str (rr)]  ,  ’ul5*u2-ul4*u3+ul3*u4-ul2*u5+ul0*u2 5 ) ; 

mem  id  =  gs  new( [’Pop’ ,num2str (rr)] ,  ’ul5*u3^-ul4*u4-ul3*u5+ull*u6+ul0*u4, ) ; 

mem  id  =  gs  new( [’Pop’ ,num2str (rr)] , ’ -ul5*u2-ul4*u3+ul3*u4-ul2*u5-ul0*u6 J ) ; 

mem  id  =  gs  new( [ ’Pop’ ,num2str (rr)] ,  ’ul5*u2+ul4*u3+ul3*u4+ul2*u5-ull*u6 ’ ) ; 

end 


This  coding  technique  builds  an  arbitrary  number,  NumPops,  of  populations,  each  “seeded” 
with  a  set  of  chromosomes  that  include  all  gene  names  and  the  arithmetic  operations  of 
addition,  subtraction,  multiplication,  and  division. 


Each  set  of  chromosomes  (one  from  each  population)  must  be  assigned  a  fitness  value. 
The  following  nest  of  loops  accomplishes  this  for  the  arbitrary  number  of  populations,  over 
an  arbitrary  number  of  receivers,  mimrec. 


mem  count  =  gs  popsize ( }Popl ’ ) ; 
for  id  =  l:mem  count 
for  rr  =  1: NumPops 

chr=strvcat (chr ,  [gs  get([,PopJ  num2str (rr)] ,  id)]); 

end 

for  rr  =  1 :  numrec 

recloc(rr, : )= [eval (chr (rr , : ))  eval(chr(rr+numrec, : )) 

eval (chr (2*numrec+l , :))  eval (chr (2*numrec+2 , :))] ; 

end 

for  z=l: numrec 

constraint (z, : )= [abs (recloc (z , 1) )  abs (recloc (z, 2))] ; 

end 

if  constraint<=[ones(numrec,l) , ones (numrec, 1)] 

fitness  =  Jfunc (recloc, xi,et a, Dlam, arc, steps, a, b,iorb,RAorb,nuorb,RAtarg,dectarg) ; 
for  rr  =  1: NumPops 

mem  id  =  gs  set  f it ( [’Pop* ,num2str (rr)] , id, fitness) ; 

end 

else 

for  rr  =  1: NumPops 

mem  del  =  gs  del(  [’Pop’ ,num2str (rr)] ,id) ; 

end 

end 

end 


The  gs  get  (...)  function  extracts  individual  chromosome  strings  from  the  population  for 
each  degree  of  freedom.  Each  is  then  converted  by  the  eval  (...)  function  into  numerical 
values.  Evaluation  of  the  Jfunc  (.  .  .)  function  is  the  main  outside  code  path.  It  calls  up 
the  master  function  that  calls  the  orbital  propagation,  calls  the  aperture  synthesis,  calls 
the  comparison  aperture  generation,  and  performs  the  aperture  comparison.  All  of  these 
calls  return  a  single  numerical  value  as  described  in  chapter  IV. 
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E.2  Performing  the  Genetic  Search 


Each  genetic  search  loop  begins  with  a  check  on  population  size.  To  ensure  a  di¬ 
verse  chromosome  population,  the  maximum  population  is  maintained  at  500.  Once  the 
population  exceeds  this,  the  gs  sel  hifit(.  .  .)  function  selects  the  character  string  of 
the  chromosome  with  the  largest  fitness  value  for  deletion  .  (There  is  a  related  function 
recommended  by  the  generic  code  writer,  gs  selr  hifit(.  .  .),  which  selects  chromo¬ 
somes  randomly  but  weighted  toward  those  with  higher  fitness  values.  Unfortunately,  the 
high  initial  “infant  mortality”  inherent  in  the  searches  made  it  nearly  inevitable  that  just 
weighting  the  selection  eventually  allowed  the  deletion  of  the  current  best  chromosome.) 


mem  count  =  gs  popsize ( ’Popl 5 ) ; 
if  mem  count  >  MaxPopulation 

mem  ids  =  gs  sel  hif it ( ’Popl ’ ) ; 
for  rr  =  l:NumPops 

gs  del(  [’Pop’ ,num2str(rr)] ,mem  ids(l)); 

end 

if  mem  count- 1  >  500 
for  rr=l:NumPops 

gs  del( [’Pop5 ,num2str (rr)] ,mem  ids(2)); 

end 

end 

end 


The  function  gs  del  (...)  actually  removes  the  chromosome  strings  from  the  populations. 


With  the  population  trimmed,  chromosomes  must  be  selected  for  the  genetic  opera¬ 
tions.  A  call  to  the  function  gs  selr  lofit(.  .  .)  selects  chromosomes  from  each  popu¬ 
lation  randomly,  but  with  the  selection  weighted  towards  those  with  lower  fitness  values. 
The  functions  gs  xover  (...)  and  gs  mut  (...)  produce  character  strings  by  concate¬ 
nating  string  portions  from  the  selected  chromosomes  to  produce  new  chromosomes.  The 
parent  chromosomes 


parentl  =  ul5/ul3+ul4*ul2  +ul3/ull+ul2/ul0+ull*ul 
parent2  =  u6*ul0+ul5+ul4  +u!3+ul2+ull 


yield  the  crossover  offspring  chromosomes 


offspringl  =  ul5/ul3+ul4*ul2  +ul3+ul2+ull 
offspring2  =  u6*ul0+ul5+ul4  +ul3/ull+ul2/ul0+ull*ul 


or  the  mutation  offspring  chromosomes 


offspringl  =  ul5/ul3+ul4*ul+ul3/ull+ul2/ul0+ull*ul2 
offspring2  =  u6*ull+ul5+ul4+ul3+ul2+ul0 


This  process  corresponds  to  step  7,  the  new  chromosome  generation  step,  of  table  6.2. 
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E.3  Evaluation  of  New  Chromosomes 

The  crossover  and  mutation  functions  generate  character  strings,  not  numerical  val¬ 
ues.  Therefore  it  is  necessary  to  again  extract  the  strings,  convert  them  to  numerical 
values,  check  them  against  the  constraints,  and  assign  fitness  values  to  them. 


for  off  =  1 : length(off  ids) 
chr=  []  ; 

for  rr  =  l:NumPops 

chr=strvcat (chr , [gs  get ( [’Pop’ ,imm2str(rr)] ,off  ids(off))]); 

end 

for  rr  =  1 :  numrec 

recloc (rr, : )= [eval (chr (rr , : ))  eval(chr(rr+numrec, : )) 

eval(chr (2*numrec+l , :))  eval (chr (2*numrec+2 , :))] ; 

end 

for  z=l: numrec 

constraint (z, : )= [abs (recloc (z , 1) )  abs(recloc(z,2))] ; 

end 

if  constraint<=[ones(numrec,l) , ones (numrec, 1)] 

fitness  =  Jfunc (recloc, xi,et a, Dlam, arc, steps, a, b,iorb,RAorb,nuorb,RAtarg,dectarg) ; 
for  rr  =  l:NumPops 

mem  id  =  gs  set  f it ( [’Pop’ ,num2str (rr)] ,off  ids (of f ) ,f itness) ; 

end 

else 

for  rr  =  l:NumPops 

mem  del  =  gs  del( [’Pop’ ,num2str (rr)] ,off  ids(off)); 

end 

end 

end 


E.4  Search  Termination 

Most  genetic  searches  execute  for  a  proscribed  number  of  generations.  The  extended 
computations  necessary  and  high  infant  mortality  in  these  searches  suggest  a  large  value  for 
NumGenerations,  but  even  that  is  insufficient  to  ensure  convergence.  Wrapping  a  while 
loop  around  the  entire  genetic  search  routine  allows  for  user  extension  of  the  search  and 
preserves  the  current  populations. 

UserContinue=’Y’ ; 

ManLoop=0; 

while  UserContinue==’ Y’ 

Genetic  Search  Loop 

UserContinue=input ( [’Want  another’ ,num2str (NumGenerations) , ’  generations?’] , ’s’) ; 
ManLoop=ManLoop+l ; 

end 

A  screen  output  (only  every  50  generations,  to  reduce  system  slow-down)  of  a  history  of 
the  best  current  fitness  value  supports  a  user  decision  of  whether  or  not  to  continue  the 
search  for  another  set  of  NumGenerations  or  allow  the  search  to  terminate. 
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Appendix  F.  Extended  Results 

This  appendix  contains  the  results  of  additional  solutions  beyond  those  directly  cited  in 
the  main  text. 


F.l  Other  General  Planar  Formations 


To  generate  the  monotonically  decreasing  plot  of  figure  6.23,  additional  optimizations 
were  performed.  Each  assumed  perfect  viewing  geometry  to  the  source. 


Convergence 
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Agg  Map  (pi/2)  2  Receivers 


Figure  F.l  A  three-receiver  formation  has  four  degrees  of  freedom,  and  a  yields  a  sur¬ 
prisingly  close  approximation  to  the  reference  aperture.  This  error  was  not 
included  in  figure  6.23  since  it  dominated  the  ordinate  axis. 
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Convergence 


Propagated  Best  Positions 


Figure  F.2  With  one  receiver  frozen  at  the  center,  a  planar  four-receiver  formation  has 
six  degrees  of  freedom. 

F.2  Other  Planar  Formations  for  Inertial  Targets 

The  inertial  targeting  algorithms  of  chapter  VII  were  further  exercised  for  larger 
formations.  Using  the  same  four  stellar  targets  as  the  five  and  six  receiver  formations 
described  in  section  7.5,  the  formations  were  expanded  to  7,  9,  and  11  receivers. 
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F.2.1  Seven  Receiver  Formations.  Expanding  the  formation  further  to  seven 
receivers  (six  movable)  yields  closer  fidelity  mimicry  of  the  filled  aperture. 
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Figure  F.3  Optimally  placing  six  of  seven  receivers  to  image  Andromeda  yields  these 
loci. 
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Figure  F.4  Reoptimizing  the  same  seven  receivers  for  Betelgeuse  yields  these  loci. 
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Figure  F.5  Seven  optimally  placed  receivers  for  viewing  Rigel  yields  these  loci. 

xIO’4  Convergence 


Loci  Viewed  from  Regulus 


Distribution  of  21  Initial  Baselines 
90  50 


Ref.  Aperture  viewed  from  Regulus 


u  (wavelengths) 
Agg  Map  (pi/2)  6  Receivers 


0  0  02 
-0.02  n  02  0 

-0.04  -0.04  . ,  ,x 

xi  (rad) 


eta  (rad) 


-0.04  -0.04" 


xi  (rad) 


Figure  F.6  Regulus’s  optimum  placement  of  six  receivers  in  the  seven  receiver  formation. 
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F.2.2  Nine  Receiver  Formations.  Expanding  the  formation  further  to  nine  re¬ 
ceivers  (eight  movable)  shows  marginal  improvements  in  synthesizing  the  filled  aperture. 
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Figure  F.7  Optimally  placing  eight  of  nine  receivers  to  image  Andromeda  yields  these 
loci. 
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Figure  F.8  Reoptimizing  the  same  nine  receivers  for  Betelgeuse  yields  these  loci. 
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Figure  F.9  Nine  optimally  placed  receivers  for  viewing  Rigel  yields  these  loci. 
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Figure  F.10 


Regulus’s  optimum  placement  of  eight  receivers  in  the  nine  receiver  forma¬ 
tion. 
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F.2.3  Eleven  Receiver  Formations.  Expanding  the  formation  further  to  eleven 


receivers  (ten  movable)  shows  marginal  improvements  in  mimicking  the  filled  aperture. 
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Figure  F.ll  Andromeda  with  eleven  (ten  movable)  receivers  has  22  degrees  of  freedom. 
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Figure  F.12 


Reoptimizing  the  formation  with  eleven  receivers  for  Betelgeuse  yields  these 
loci. 
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Figure  F.13  Eleven  optimally  placed  receivers  for  viewing  Rigel  yields  these  loci. 
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Figure  F.14  Regulus’s  optimum  placement  of  ten  receivers  in  the  eleven  receiver  forma¬ 
tion. 
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